Supporting Code

for “Animal-vectored nutrient flows over resource gradients influence the nature of local and meta-ecosystem functioning”

Introduction

This notebook contains the code used to run numerical analyses of a meta-ecosystem model in which consumers move across ecosystem borders in different ways with respect to resource gradients.

In section Model Numerical Analyses, we perform numerical analyses of our model of consumer movement in a two-patch meta-ecosystem, using the set of equilibrium equations derived in Mathematica (see Mathematica notebook in this same repository; Wolfram Research Inc. (2020)). In section Changes to biomass distribution in local and meta-ecosystem, we investigate how biomass distribution changes following different types of consumer movement and how this may affect top-down influences of consumer movement on local and meta-ecosystem functions. Finally, in section Sensitivity Analyses, we perform a suite of sensitivity analyses on our results.

In this introduction, we briefly define the contrasting consumer movement types we investigate with this model (Defining movement types), discuss separation of scales in time and space (Separation of Spatial and Temporal scales), describe the state variables and parameters in the model (State variables and parameters of the model), and describe the analyses setup (Analyses setup).

Defining movement types

Here, we focus on consumers moving over gradients of inorganic nutrients availability between two local ecosystems. Consumer movement—either gradient neutral, along-, or against-gradient—connects these two local ecosystems in a meta-ecosystem (Loreau, Mouquet, and Holt 2003). We define consumer movement as gradient neutral when it happens among ecosystems that have equal nutrient availability. Along-gradient movement takes place when consumers move following the gradient direction. Conversely, against-gradient consumer movement runs counter to the direction of the gradient. Consistent with recent reviews (Gounand et al. 2018), in the following sections we also use diffusive and non-diffusive to refer to along- and against-gradient consumer movement, respectively. In our model, we create resource gradients by varying the relative values of the inorganic nutrient input rates in the two ecosystems comprising our meta-ecosystem (parameter I_i; see Table 1 in the main text or below for a list of model parameters and state variables).

Separation of spatial and temporal scales

While organisms can at times cross into an ecosystem that is equally suitable for them immediately after leaving their original one, this is not often the case. Rather, organisms spend significant amounts of time searching for a new suitable ecosystem, traversing the so-called “unsuitable” matrix (Weisser and Hassell 1996). Models of organismal movement in meta-ecosystems should account for this time spent outside ecosystem borders, as it may lead to loss of biomass and nutrients from the system (Gounand et al. 2018).

We model consumer movement in our meta-ecosystem model by combining a dispersers’ pool (Weisser and Hassell 1996) with the time scales separation (TSS) mathematical technique (Otto and Day 2011). Briefly, adding a dispersers’ pool to a classic meta-ecosystem model means that consumer do not instantly travel from donor to recipient ecosystem. Instead, they move from the donor to the dispersers’ pool and then from there to the recipient (see Figure 1c in the main text for a visual representation). This allows us to separate ecosystem processes that happen locally from those that happen regionally.

Some of these ecosystem processes, like primary productivity in terrestrial ecosystems, take place over long time frames (slow processes). Others, like consumer movement, happen over short time frames (fast processes). TSS explicitly accounts for these differences by allowing fast processes to reach a quasi-equilibrium while holding slow ones constant. The quasi-equilibrium is then plugged into the equation(s) for the slow processes to investigate how these processes evolve as the fast ones stay constant (Otto and Day 2011).

The combined use of the dispersers’ pool and TSS allows us to separate local and regional dynamics in our system and to investigate the effects of consumer movement at either spatial extent at the same time. See section Model Derivation in the main text for further details on our novel approach to model consumer movement in meta-ecosystems.

State variables and parameters of the model

The table below lists the model’s state variables and parameters, and is analogous to Table 1 in the main text. For a visual representation of the model, please see Figure 1c in the main text.

Variable	Description	Units	Range
Ni	Nutrients stocks in patch i	g	>0
Pi	Producers stocks in patch i	g	>0
Ci	Consumers stocks in patch i	g	>0
Q	Dispersers’ Pool	g	>0

Parameter	Description	Units	Range
Ii	Inorganic nutrient input rate into patch i	g*t-1	[0,20]
\(l\)	Inorganic nutrient output rate	t-1	[0,10]
ui	Producer uptake rate in patch i	(g*t)-1	[0,10]
ai	Consumer attack rate in patch i	(g*t)-1	[0,10]
\(\epsilon\)i	Consumer assimilation efficiency in patch i	unitless	[0,1]
hi	Biomass loss rate from Pi	t-1	[0,10]
di	Biomass loss rate from Ci	t-1	[0,10]
g	Movement rate from C1 to Q	t-1	[0,10]
m	Movement rate from Q to C2	t-1	[0,10]
c	Biomass loss rate from Q	t-1	[0,10]

We derived the model’s single feasible equilibrium in Mathematica (Wolfram Research Inc. (2020); see Model Derivation in the main text for further details).

Analyses setup

This notebook uses the following packages:

Show code

# The following packages are used throughout this document.  Before
# loading each packages, we provide aa brief description of what it
# is used for

# project dependency (packages, mostly) management
library(renv)
# data wrangling, manipulation, plotting
library(tidyverse)
# additional plotting functions
library(plotrix)
# use latin hypercube to create toy datasets
library(lhs)
# building tables in a tidyverse framework
library(gt)
# scaling functions for ggplot2
library(scales)
# color-blind palettes for ggplot2
library(khroma)
# extra tabbing functionality for rmarkdown documents
library(xaringanExtra)
# functions to time compilation of this document
library(tictoc)
# color palettes manipulation
library(RColorBrewer)
# additional color palettes inspired by MET artworks
library(MetBrewer)
# additional data wrangling
library(broom)
# print local directory structure
library(fs)
# additional ggplot2 themes
library(cowplot)
# use randomForest algorithms to perform sensitivity analyses on
# results
library(randomForest)
# extra themes for ggplot2
library(ggpubr)
# extra geoms for ggplot2
library(ggpol)
# easily create and layout multi-panel plots
library(patchwork)
# additional ggplot2 functionality
library(ggh4x)
# extra themes for ggplot2
library(ggpubr)
# extra geoms for ggplot2
library(ggpol)
# easily create and layout multi-panel plots
library(patchwork)

# Packages that need to be installed but do not require loading

# library(distill) # notebook templates library(rmarkdown) # build
# interactive documents in R

# set the theme for all ggplot2 objects
theme_set(theme_classic())

In the interest of reproducibility, we use package renv to keep track of this notebook’s packages.

Show code

renv::snapshot()

* Lockfile written to '~/PhD/Chapter_3/Code/renv.lock'.

We compiled this notebook on a machine using the following version of R, operative system, and necessary packages:

Show code

# this needs three colons to work between package name and function call
utils:::print.sessionInfo(sessionInfo()[-8])

R version 4.2.2 (2022-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.6.3

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods  
[7] base     

other attached packages:
 [1] ggh4x_0.2.3          patchwork_1.1.2      ggpol_0.0.7         
 [4] ggpubr_0.5.0         randomForest_4.7-1.1 cowplot_1.1.1       
 [7] fs_1.6.0             broom_1.0.3          MetBrewer_0.2.0     
[10] RColorBrewer_1.1-3   tictoc_1.1           xaringanExtra_0.7.0 
[13] khroma_1.9.0         scales_1.2.1         gt_0.8.0            
[16] lhs_1.1.6            plotrix_3.8-2        forcats_1.0.0       
[19] stringr_1.5.0        dplyr_1.1.0          purrr_1.0.1         
[22] readr_2.1.3          tidyr_1.3.0          tibble_3.1.8        
[25] ggplot2_3.4.0        tidyverse_1.3.2      renv_0.16.0         

Finally, this notebook assumes the following folder structure:

Show code

dir_tree(path = "../", recurse = 0, all = FALSE, type = "directory")

../
├── Code
├── Data
├── LitRev
├── Manuscripts
├── Modeling_Ch3
├── Results
└── Submission

This structure allows us to take advantage of relative paths when loading data or saving figures. Note that ../ stands for the root folder. If you wish to recompile this notebook, please organize your folder structure accordingly.

Model numerical analyses

Here, we perform the analyses of the meta-ecosystem model described in section Numerical Analyses of the main text. If you are interested in the code used to produce the main text’s figures and tables, this can be found in sections Visual Deliverables—Ecosystem functions and Summary Tables—Ecosystem functions below.

We conduct three rounds of analyses, one each for the three types of consumer movement we are interested in: gradient-neutral, along-, and against-gradient. Each round of analysis comprises iterating the model over n = 10 000 parameter sets, running stability analyses, and graphing the results. The analyses will follow these steps:

1. Fit randomly-drawn parameter sets to the model’s equilibrium and solve for state variable values
2. Calculate meta-ecosystem function values
3. Flag parameter sets that produce equilibria < 0, hence without biological meaning
4. Run stability analyses on each parameter set to identify those producing unstable results
5. Exclude parameter sets that are unstable and/or have no biological meaning from further analyses
6. Explore results graphically

Generating parameter sets

This notebooks defaults to loading the two pre-generated datasets provided alongside it, in folder ../Data. DATA1 includes parameters scaled [0–10], whereas DATA2 includes parameters scaled [0–1]. These were generated with a Latin hypercube sampling design, using function lhs(), and then used to produce the results and figures presented in the paper. To generate new randomly drawn parameter values, set eval=FALSE in the data-load code chunk below and eval=TRUE in the following one, labelled disp-poolsetParams.

Show code

DATA1 <- read.csv(file = "../Data/DATA1.csv", header = TRUE)

DATA2 <- read.csv(file = "../Data/DATA2.csv", header = TRUE)

In the code chunk below, we show how to generate two sets of randomly-drawn parameter values to use in all scenarios. We use function lhs() to do so. One set includes parameters values scaled 0–10 (DATA1), while the other one includes parameter values scaled 0–1 (DATA2).

Show code

# Below we generate random values for those parameters that will vary
# during the simulations: attack rates, efficiency, death rates, and
# movement rates

DATA1 <- NULL
# Run 2000 iterations of 100 numbers (in 100 bins between 0 to 1)
for (i in seq(1, 100, 1)) {

    # Use lhs function in R. 100 bins, 11 iterations - one for each
    # parameter that is scaled between 0 and 10.
    loop_1 <- randomLHS(100, 11)

    # Use rbind to stack each iteration on top of the other to create
    # a one column list of 10000 samples with equal distribution in
    # 100 bins.  Multiplied by 10 to get range of data from 0 to 10
    DATA1 <- rbind(DATA1, data.frame(i = i, Result = loop_1 * 10))

    # print the i value so you can track progress print(i) Close the
    # loop
}

# This is for 2 parameters that is scaled 0 to 1.
DATA2 <- NULL
# Run 2000 iterations of 100 numbers (in 100 bins between 0 to 1)
for (i in seq(1, 100, 1)) {

    # Use lhs function in R. 100 bins, 2 iterations.
    loop_2 <- randomLHS(100, 2)

    # Use rbind to stack each iteration on top of the other to create
    # a one column list of 200000 samples with equal distribution in
    # 100 bins.  Multiplied by 10 to get range of data from 0 to 10
    DATA2 <- rbind(DATA2, data.frame(i = i, Result = loop_2))

    # print the i value so you can track progress print(i) Close the
    # loop
}

# Save these two datasets for later re-use write.csv(DATA1, file =
# '../Data/DATA1.csv', row.names = FALSE) write.csv(DATA2, file =
# '../Data/DATA2.csv', row.names = FALSE)

Control scenario: gradient-neutral consumer movement

We begin our numerical analyses of our model with the control scenario, in which consumers move between two ecosystems that have equal environmental fertility. Consumer movement in this scenario is gradient-neutral, as the environment is homogeneous in terms of resource availability. The results of this first set of model iterations produce the values used at the denominator of the response ratio (LRR) used to capture changes in local and meta-ecosystem function values (see section [Visual Deliverables] below). Below, we describe each step of the analyses in detail: later analyses will not describe these steps again.

Step 1: plug the parameter sets into the model

First, we allocate an empty dataframe PRED to contain the results of each iteration of the model. This dataframe will contain the parameter values and the resulting local and meta-ecosystem function values.

Show code

# Allocate an empty data frame to save the data in, then allocate
# some rows and columns. We will replace the current numbers with
# results below Data we are calculating are (bio)mass for soil
# nutrients (N1, N2), producers (P1, P2), consumers (C1, C2),
# disperses (Q)

PRED <- NULL

PRED <- rbind(PRED, data.frame(TIME = seq(0, 100, 1), I1 = seq(0, 100,
    1), I2 = seq(0, 100, 1), l = seq(0, 100, 1), u1 = seq(0, 100, 1), u2 = seq(0,
    100, 1), a1 = seq(0, 100, 1), a2 = seq(0, 100, 1), h1 = seq(0, 100,
    1), h2 = seq(0, 100, 1), d1 = seq(0, 100, 1), d2 = seq(0, 100, 1),
    g = seq(0, 100, 1), m = seq(0, 100, 1), c = seq(0, 100, 1), e1 = seq(0,
        100, 1), e2 = seq(0, 100, 1), N1 = seq(0, 100, 1), P1 = seq(0,
        100, 1), C1 = seq(0, 100, 1), N2 = seq(0, 100, 1), P2 = seq(0,
        100, 1), C2 = seq(0, 100, 1), Q = seq(0, 100, 1), FLUX_P1 = seq(0,
        100, 1), FLUX_C1 = seq(0, 100, 1), FLUX_P2 = seq(0, 100, 1), FLUX_C2 = seq(0,
        100, 1), FLUX_Eco_1 = seq(0, 100, 1), FLUX_Eco_2 = seq(0, 100,
        1), FLUX_Q = seq(0, 100, 1), PROD_P1 = seq(0, 100, 1), PROD_C1 = seq(0,
        100, 1), PROD_P2 = seq(0, 100, 1), PROD_C2 = seq(0, 100, 1)))

Now we run the model. We first set the values for both nutrient input rates in ecosystem 1 (I1 in the code chunks) and ecosystem 2 (I2) at 10. Recall that this is the control scenario with equal environmental fertility. We also set the value of the inorganic nutrient leaching rate (l) at 0.1.

Note that we refer to I1 and I2 as either the nutrient input rates or the environmental fertility of the two ecosystems

The code chunk below runs the model using a for loop that, over 10 000 time steps, iterates the model’s single feasible equilibrium to calculate:

• stocks of Nutrients (N1, N2), Producers (P1, P2), and Consumers (C1, C2)
• the Dispersers’ Pool (Q) stock
• the recycling flux for each biotic compartment (FLUX_X, where X is the compartment)
• the productivity for each biotic compartment (PROD_X, where X is the trophic compartment)

in both local ecosystems. At each time step, we draw a random value for each parameter from either the DATA1 or DATA2 dataset—depending on whether the parameter is scaled 0–10 or 0–1, respectively. Then, we plug these values into the model’s equilibria and ecosystem function formulas derived in Mathematica (Wolfram Research Inc. 2020), and finally store both the randomly drawn parameter values and ecosystem function values in PRED.

Note that feasibility conditions for this model cannot be resolved analytically. Hence, we do not check that our random parameters values satisfy them. Rather, we run stability analyses of each parameter set below.

Show code

# Here we set values for parameters that we will not vary during the simulations
# Nutrients input rate in each patch
I1 = 10
I2 = 10
# Patch leaching rate/soil loss rate
l = 0.1

# We want to start simulations at 0
i1 = 0

for (i in seq(0,9999,1)){
  # go to the next row
  i1 = i1 + 1

# Let's sample the parameter values from DATA. We will sample these i times as 
# defined by the for loop above

  # producers uptake rates
  u1 = DATA1[i1,2] # in ecosystem 1
  u2 = DATA1[i1,3] # in ecosystem 2
  # consumers attack rates
  a1 = DATA1[i1,4] # in ecosystem 1
  a2 = DATA1[i1,5] # in ecosystem 2
  # death rates (= recycling rate)
  h1 = DATA1[i1,6] # death (=recycling) rate for producers in ecosystems 1
  h2 = DATA1[i1,7] # death (=recycling) rate for producers in ecosystems 2
  d1 = DATA1[i1,8] # death (=recycling) rate for consumers in ecosystem 1
  d2 = DATA1[i1,9] # death (=recycling) rate for consumers in ecosystem 2
  c = DATA1[i1,10] # death rate in the disperser's pool Q
  # movement rate to the Dispersers' pool Q
  g = DATA1[i1,11]
  # movement rate from the Dispersers' pool Q
  m = DATA1[i1,12]
  # assimilation efficiency of consumers 
  e1 = DATA2[i1,2] # in ecosystem 1
  e2 = DATA2[i1,3] # in ecosystem 2
  
  PRED[i1,] <- c(i1, I1, I2, l, u1, u2, a1, a2, h1, h2, d1, d2, g, m, c, e1, e2,
                 # Nutrient stock in ecosystem 1
                 ((d1 - d1*e1 + g)*h1 + a1*e1*I1)/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 # Primary Producers biomass in ecosystem 1
                 (d1 + g)/(a1*e1),
                 # Consumers biomass in ecosystem 1
                 (e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 # Nutrient stock in ecosystem 2
                 (-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)),
                 # Primary Producers biomass in ecosystem 2
                 (l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
                 # Consumer biomass in ecosystem 2
                 ((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2))/(a2*e2*l - d2*(-1 + e2)*u2),
                 # Consumer biomass in the Dispersers' Pool (Q)
                 (e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
                 # Nutrient Flux for Primary Producers in Ecosystem 1
                 ((d1 + g)*h1)/(a1*e1),
                 # Nutrient Flux for Consumers in Ecosystem 1
                 (d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 # Nutrient Flux for Primary Producers in Ecosystem 2
                 (h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
                 # Nutrient Flux for Consumers in Ecosystem 2
                 (d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
                 # Nutrient Flux for Ecosystem 1
                 ((d1 + g)*h1)/(a1*e1) + (d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 # Nutrient Flux for Ecosystem 2
                 (h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2) + (d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
                 # Nutrient Flux in the Dispersers' Pool (Q)
                 (c*e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
                 # Primary Productivity in Ecosystem 1
                 ((d1 + g)*((d1 - d1*e1 + g)*h1 + a1*e1*I1)*u1)/(a1*e1*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
                 # Consumer (=Secondary) Productivity in Ecosystem 1
                 (e1*(d1 + g)*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 # Primary Productivity in Ecosystem 2
                 ((-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))*u2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)*(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2)),
                 # Consumer (=Secondary) Productivity in Ecosystem 2
                 (e2*l*(-a1*d2*e1*h2*l*(c + m) + d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(h1*l - I1*u1)) + d2*e2*(a1*e1*I2*l*(c + m) + (d1 + g)*I2*(c + m)*u1 - e1*(d1*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2))
                 )
}

Step 2: calculate meta-ecosystem functions values

To investigate any influence of consumer movement on nutrient cycling at regional spatial scales, we calculate local and meta-ecosystem total biomass (Ntot,Ptot, Ctot), meta-ecosystem nutrient flux (MetaEcoFlux) and trophic compartment nutrient flux (FLUX_Ptot, FLUX_Ctot) and productivity (PROD_Ptot, PROD_Ctot). Furthermore, we manually calculate each ecosystem’s own total recycling flux as the sum of the flux of its 2 biotic compartments, as a check on the Mathematica-derived formulas used above.

Show code

# calculate meta-ecosystem biomass stocks for each trophic
# compartment, a.k.a., meta-ecosystem biomass
PRED$Ntot <- PRED$N1 + PRED$N2
PRED$Ptot <- PRED$P1 + PRED$P2
PRED$Ctot <- PRED$C1 + PRED$C2

# now calculate recycling flux for the local ecosystem manually to
# double check the formulae used above are correct
PRED$FLUX_Eco1_check <- PRED$FLUX_P1 + PRED$FLUX_C1
PRED$FLUX_Eco2_check <- PRED$FLUX_P2 + PRED$FLUX_C2

# and the meta-ecosystem recycling flux
PRED$FLUX_Ptot <- PRED$FLUX_P1 + PRED$FLUX_P2
PRED$FLUX_Ctot <- PRED$FLUX_C1 + PRED$FLUX_C2
PRED$MetaEcoFlux <- PRED$FLUX_Eco_1 + PRED$FLUX_Eco_2 - PRED$FLUX_Q

# finally, calculate the production for each compartment at
# meta-ecosystem level
PRED <- dplyr::mutate(PRED, PROD_Ptot = PROD_P1 + PROD_P2, .after = "PROD_C2")
PRED <- dplyr::mutate(PRED, PROD_Ctot = PROD_C1 + PROD_C2, .after = "PROD_Ptot")

Step 3: flag parameter sets with no biological meaning

Parameter sets that produce values of the model’s state variable that are \(\leq 0\) lack biological meaning. Accordingly, in the code chunk below, we flag these parameter sets for later exclusion from the analysis.

Show code

# Equilibrium stocks less than or equal to 0 make no biological
# sense, so we will flag them in the PRED dataset.
PRED <- PRED %>%
    dplyr::mutate(., biosense = ifelse(N1 > 0 & P1 > 0 & C1 > 0 & N2 >
        0 & P2 > 0 & C2 > 0 & Q > 0, "yes", "no"), .after = MetaEcoFlux) %>%
    dplyr::mutate_at(vars(biosense), factor)

Step 4: stability analyses

We perform the stability analyses on the model by fitting the same sets of parameter values and equilibrium state variable values used above to the partial derivatives extracted from the model’s Jacobian matrix evaluated in Mathematica (Wolfram Research Inc. 2020). We evaluate stability of each parameter set by checking whether the real part of the leading eigenvalue of the Jacobian is positive or negative (Otto and Day 2011). If the real part of the Jacobian’s leading eigenvalue is negative, the equilibrium is stable (Otto and Day 2011).

Show code

# create an empty dataframe
MathStab <- NULL

for (i in 1:nrow(PRED)) {

    # fetch parameter values from the equilibrium simulations df
    I1 = PRED$I1[i]
    I2 = PRED$I2[i]
    l = PRED$l[i]
    u1 = PRED$u1[i]
    u2 = PRED$u2[i]
    a1 = PRED$a1[i]
    a2 = PRED$a2[i]
    h1 = PRED$h1[i]
    h2 = PRED$h2[i]
    d1 = PRED$d1[i]
    d2 = PRED$d2[i]
    g = PRED$g[i]
    m = PRED$m[i]
    c = PRED$c[i]
    e1 = PRED$e1[i]
    e2 = PRED$e2[i]

    # Put in the solutions to the equilibrium values here:
    dN1 = PRED$N1[i]
    dP1 = PRED$P1[i]
    dC1 = PRED$C1[i]
    dN2 = PRED$N2[i]
    dP2 = PRED$P2[i]
    dC2 = PRED$C2[i]

    # create the Jacobian from Mathematica NOTE: in transposing from
    # Mathematica, the Jacobian is here assembled by row
    Jacob <- rbind(c(-l - dP1 * u1, h1 - dN1 * u1, d1, 0, 0, 0), c(dP1 *
        u1, -a1 * dC1 - h1 + dN1 * u1, -a1 * dP1, 0, 0, 0), c(0, a1 * dC1 *
        e1, -d1 - g + a1 * e1 * dP1, 0, 0, 0), c(0, 0, 0, -l - dP2 * u2,
        h2 - dN2 * u2, d2), c(0, 0, 0, dP2 * u2, -a2 * dC2 - h2 + dN2 *
        u2, -a2 * dP2), c(0, 0, (g * m)/(c + m), 0, a2 * dC2 * e2, -d2 +
        a2 * e2 * dP2))

    # Check to see if the real part of the leading eigenvalue is less
    # than 0 or not - if it is than the system is stable.
    if (max(Re(base::eigen(Jacob)$values)) < 0) {
        stable <- "stable"
    } else {
        stable <- "unstable"
    }

    MathStab <- rbind(MathStab, data.frame(TIME = i, I1, I2, l, u1, u2,
        a1, a2, h1, h2, d1, d2, g, m, c, e1, e2, dN1, dP1, dC1, dN2, dP2,
        dC2, EiV1 = Re(eigen(Jacob)$values[1]), EiV2 = Re(eigen(Jacob)$values[2]),
        EiV3 = Re(eigen(Jacob)$values[3]), EiV4 = Re(eigen(Jacob)$values[4]),
        EiV5 = Re(eigen(Jacob)$values[5]), EiV6 = Re(eigen(Jacob)$values[6]),
        maxEv = max(Re(base::eigen(Jacob)$values)), stable = stable, biosense = PRED$biosense[i]))
}

MathStab$stable <- as.factor(MathStab$stable)

# separate unstable equilibria to work with later
MathStabUS <- subset(MathStab, MathStab$stable == "unstable")

According to the stability analyses run above, about 1.55% of the parameter sets produced unstable results (i.e., 155 out of 10000 iterations).

Step 5: exclude unstable parameter sets from further analyses

Now, let’s check if the parameter sets that produce unstable equilibria are also those that produce equilibria lacking biological meaning—i.e., those where state variables values at equilibrium are \(\leq 0\).

Show code

biononsense <- subset(PRED, PRED$biosense == "no")

identical(as.numeric(MathStabUS[, "TIME"]), biononsense[, "TIME"])

[1] TRUE

It appears that they are. Hence, let’s add the stable/unstable information to PRED.

Show code

PRED <- left_join(PRED, select(MathStab, !c(I1:maxEv)), by = c("TIME",
    "biosense"))

Now, let’s exclude these unstable parameter sets from the analyses below and from our results.

We are also going to sample the dataset PRED, that contains the results of our model’s iterations, for 1000 random iterations results. We will store these in PRED_1k and use them in Figures 1, 2, 3, 4, 5 below to produce some preliminary graphs of our results.

Show code

PREDpos <- filter(PRED, PRED$biosense == "yes" & PRED$stable == "stable")

# sample PREDpos to only use 1000 random simulation results
PRED_sample1000 <- PREDpos[sample(nrow(PREDpos), size = 1000), ]

PRED_1k <- droplevels(PRED_sample1000)

Before proceeding further, let’s save the PRED object to disk, for ease of use in the future.

Show code

write.csv(PRED, file = "../Results/PRED.csv", row.names = FALSE)

Step 6: explore results graphically

We now investigate how the presence of unidirectional, gradient-neutral consumer movement influences the meta-ecosystem’s behavior using box-plots produced using PRED_1k. Note that in the figures below, “Donor” always indicates ecosystem 1 and “Recipient” always indicates ecosystem 2.

Biomass and nutrient stock accumulation

Show code

# pivot the dataset to longer format for boxplot plotting
PRED_biomass_long <- select(PRED_1k, N1:Q, Ntot:Ctot) %>%
    pivot_longer(., names_to = "Compartment", values_to = "Stock", cols = c(N1,
        P1, C1, N2, P2, C2, Q, Ntot, Ptot, Ctot)) %>%
    dplyr::mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
        "C1" | Compartment == "P1", "Donor", if_else(Compartment == "N2" |
        Compartment == "P2" | Compartment == "C2", "Recipient", if_else(Compartment ==
        "Q", "Dispersers Pool", "Meta-ecosystem"))))

PRED_biomass_long$Compartment <- as_factor(PRED_biomass_long$Compartment)
PRED_biomass_long$Ecosystem <- as_factor(PRED_biomass_long$Ecosystem)

comparts_medians <- PRED_biomass_long %>%
    dplyr::group_by(Compartment, Ecosystem) %>%
    dplyr::summarise(., median = median(Stock), .groups = "keep")

ggplot(data = PRED_biomass_long, aes(x = Compartment, y = Stock, col = Ecosystem)) +
    geom_boxjitter(alpha = 0.25, outlier.intersect = TRUE, outlier.shape = 25,
        jitter.shape = 21, outlier.stroke = 1, jitter.stroke = 1) + stat_summary(geom = "label",
    fun = quantile, fun.args = list(probs = 0.5), aes(label = round(after_stat(y),
        2)), position = position_dodge(width = 0.2), show.legend = F, size = 3,
    label.padding = unit(0.15, "lines"), label.size = 0.15, fontface = "bold",
    color = "black") + facet_grid(. ~ Ecosystem, scales = "free") + scale_color_manual(values = met.brewer("Egypt",
    4)) + theme_pubr() + coord_cartesian(y = c(0, 100)) + theme(legend.position = "none")

Figure 1: Biomass stock values at equilibrium for local and meta-ecosystem scale, drawn using 1k randomly-drawn data points from the model simulations. Consumer movement from ecosystem 1 to ecosystem 2 via the dispersers’ pool results in foraging pressure of consumers being transferred from producers in ecosystem 1 to those in ecosystem 2. Primary producers in ecosystem 1 are thus released from foraging pressure while those in ecosystem 2 are depressed. In turn, this leads to differences in nutrient stocks between the two ecosystems at both the local and meta-ecosystem extents. Point-down triangles identify outliers; labels on top of the compartments’ boxplots report median values for the full set of simulations (i.e., 10k iterations). Values on the y-axis are limited between [0, 100] to zoom in on the variation in the results.

Show code

# ggsave(filename = '../Results/ModelEq_base.pdf', device = 'pdf',
# dpi = 300, height = 5)

The graph above shows a spatial trophic cascade (Knight et al. 2005). In ecosystem 2, the influx of consumers from ecosystem 1 increases the consumer biomass, in turn depressing the producers’ abundance and releasing the nutrient stock from the producers’ trophic pressure. Conversely, in ecosystem 1, consumers’ numbers fall steadily—due to the combined action of mortality and emigration towards ecosystem 2—thus releasing the producers from the trophic pressure exerted by consumers. In turn, this lets producers grow unchecked, depressing the nutrient stocks in ecosystem 1. We can see this as well by looking at the change of log₁₀ biomass values for each trophic compartment with respect to the two movement parameters: g, the rate of consumer movement from ecosystem 1 to the dispersers’ pool Q (Figure 2), and m, the rate of consumer movement from the dispersers’ pool Q to ecosystem 2 (Figure 3).

Show code

eco_levels <- c(N1 = "Nutrient", P1 = "Producers", C1 = "Consumers", N2 = "Nutrient",
    P2 = "Producers", C2 = "Consumers", Q = "Dispersers", Ntot = "Nutrient",
    Ptot = "Producers", Ctot = "Consumers")

# pivot the dataset to longer format for boxplot plotting
PRED_long <- pivot_longer(PRED_1k, names_to = "Compartment", values_to = "Stock",
    cols = c(N1, P1, C1, N2, P2, C2, Q, Ntot, Ptot, Ctot)) %>%
    dplyr::mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
        "C1" | Compartment == "P1", "Donor", if_else(Compartment == "N2" |
        Compartment == "P2" | Compartment == "C2", "Recipient", if_else(Compartment ==
        "Q", "Dispersers", "Meta-ecosystem"))))

PRED_long$Compartment <- as_factor(PRED_long$Compartment)
PRED_long$Ecosystem <- as_factor(PRED_long$Ecosystem)

ggplot(PRED_long, aes(g, log10(Stock), col = Ecosystem)) + geom_point(alpha = 0.15,
    shape = 21, size = 2, stroke = 1) + geom_smooth(method = "lm", col = "white",
    se = T, lwd = 0.4) + facet_grid(. ~ Compartment) + scale_color_manual(values = met.brewer("Egypt",
    4)) + theme_pubr() + guides(colour = guide_legend(override.aes = list(alpha = 1))) +
    theme(axis.text.x = element_text(angle = 45, vjust = 0.5, size = 8),
        text = element_text(size = 11))

Figure 2: Change in biomass in the meta-ecosystems’ trophic compartments with change in the movement rate from the donor ecosystem 1 to the dispersers’ pool (g). Each facet presents results for a trophic compartment, with color identifying local and meta-ecosystem. White lines are regression lines with 95% CI (grey shaded areas). Note the log₁₀ scale on the y-axis. Graph produced using 1000 randomly-drawn data points from the model’s analyses.

Show code

ggplot(PRED_long, aes(m, log10(Stock), col = Ecosystem)) + geom_point(alpha = 0.15,
    shape = 21, size = 2, stroke = 1) + geom_smooth(method = "lm", col = "white",
    se = T, lwd = 0.4) + facet_grid(. ~ Compartment) + scale_color_manual(values = met.brewer("Egypt",
    4)) + theme_pubr() + guides(colour = guide_legend(override.aes = list(alpha = 1))) +
    theme(axis.text.x = element_text(angle = 45, vjust = 0.5, size = 8),
        text = element_text(size = 11))

Figure 3: Change in biomass in the meta-ecostystem’s trophic compartments with change in movement rate from the disperser pool to the recipient ecosystem 2 (m). Graph produced using 1000 randomly-drawn data points from the model’s analyses. Note the log₁₀ scale on the y-axis. All other specifications as in Figure 2.

Finally, we can see a similar effect when looking at the change in biomass in the two local ecosystems with change in the mortality rate of consumers while in the dispersers’ pool (parameter c in the model; Figure 4).

Show code

ggplot(PRED_long, aes(c, log10(Stock), col = Ecosystem)) + geom_point(alpha = 0.15,
    shape = 21, size = 2, stroke = 1) + geom_smooth(method = "lm", col = "white",
    se = T, lwd = 0.4) + facet_grid(. ~ Compartment) + scale_color_manual(values = met.brewer("Egypt",
    4)) + theme_pubr() + guides(colour = guide_legend(override.aes = list(alpha = 1))) +
    theme(axis.text.x = element_text(angle = 45, vjust = 0.5, size = 8),
        text = element_text(size = 11))

Figure 4: Change in biomass in the meta-ecostystem’s trophic compartments with change in the consumer death rate while in the dispersers’ pool (c). This could represent predation risk, but also environmental hazards and stochastic events such as crossing a dangerous linear feature of the landscape—e.g., a busy highway. Graph produced using 1000 randomly-drawn data points from the model’s analyses. Note the log₁₀ scale on the y-axis. All other specifications as in Figure 2.

Nutrient flux

In the graphs below (Figure 5), we look at nutrient flux values in each ecosystem and in the meta-ecosystem for each iteration of the model.

First, we will extract the flux columns (see above) from the PRED_1k dataframe into a new one which we will call FluxPRED. We will then pivot this dataframe from wide to long format, using function pivot_longer from package tidyr.

Show code

FluxPRED <- PRED_1k %>%
    dplyr::select(TIME:c, FLUX_P1:FLUX_Eco_2, MetaEcoFlux) %>%
    tidyr::pivot_longer(., names_to = "Compartment", values_to = "Flux",
        cols = c(FLUX_P1:MetaEcoFlux)) %>%
    dplyr::mutate(., Scale = if_else(Compartment == "FLUX_P1" | Compartment ==
        "FLUX_C1" | Compartment == "FLUX_Eco_1", "Donor", if_else(Compartment ==
        "FLUX_P2" | Compartment == "FLUX_C2" | Compartment == "FLUX_Eco_2",
        "Recipient", "Meta-ecosystem"))) %>%
    dplyr::mutate(., Scale = fct_relevel(Scale, "Donor", "Recipient", "Meta-ecosystem"))

Now, we produce the box-and-whiskers plot of recycling flux in local and meta-ecosystem.

Show code

FluxPRED$Scale <- as.factor(FluxPRED$Scale)
FluxPRED$Compartment <- as.factor(FluxPRED$Compartment)

comparts_medians <- FluxPRED %>%
    dplyr::group_by(Compartment, Scale) %>%
    dplyr::summarise(., median = median(Flux), .groups = "keep")

FluxPRED %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, c("FLUX_P1",
        "FLUX_C1", "FLUX_Eco_1", "FLUX_P2", "FLUX_C2", "FLUX_Eco_2", "MetaEcoFlux"))) %>%
    dplyr::mutate(., Scale = fct_relevel(Scale, c("Donor", "Recipient",
        "Meta-ecosystem"))) %>%
    ggplot(., aes(x = Compartment, y = Flux)) + geom_boxjitter(outlier.intersect = TRUE,
    alpha = 0.25, outlier.shape = 25, jitter.shape = 21) + stat_summary(geom = "label",
    fun = quantile, fun.args = list(probs = 0.5), aes(label = round(after_stat(y),
        2)), position = position_dodge(width = 0.2), show.legend = F, size = 3,
    label.padding = unit(0.15, "lines"), label.size = 0.15, fontface = "bold",
    color = "black") + facet_grid(. ~ Scale, scales = "free") + scale_x_discrete(labels = c(FLUX_C1 = "C1",
    FLUX_P1 = "P1", FLUX_Ptot = "Producers", FLUX_C2 = "C2", FLUX_P2 = "P2",
    FLUX_Ctot = "Consumers", FLUX_Eco_1 = "Donor", FLUX_Eco_2 = "Recipient",
    MetaEcoFlux = "Meta-ecosystem")) + theme_pubr() + coord_cartesian(y = c(0,
    500)) + theme(legend.position = "none", axis.title.x = element_blank())

Figure 5: Nutrient flux at local and meta-ecosystem scales when consumers move from ecosystem 1 to 2. The two ecosystems have the same fertility—i.e., basal inorganic input values are the same. Point-down triangles identify outliers; labels on top of the compartments’ boxplots report median values for the full set of simulations. The y-axis is constrained between [0, 500] to show the variation in the results. Graph produced using 1000 randomly-drawn data points from the model’s analyses.

Experimental scenarios: along- and against-gradient consumer movement

Here, we investigate the influence of consumers moving from ecosystem 1 to ecosystem 2 when there is a gradient of nutrient availability between the two ecosystems. We will investigate two scenarios: higher environmental fertility in the donor ecosystem 1 (i.e., I1 >> I2) and higher environmental fertility in the recipient ecosystem 2 (i.e., I1 << I2). Note that the environmental leaching rate of each ecosystem—i.e., the rate at which nutrients leave local ecosystems independently of trophic processes—is invariant and equal between the two ecosystems (l = 0.1).

The analysis will follow the same steps as before: draw parameter values from our DATA1 and DATA2 data objects, simulate model behavior, run stability analyses, remove equilibria that are unstable and/or lack biological sense, and graphically inspect and interpret the results.

Along-gradient consumer movement

We begin by simulating a meta-ecosystem where the donor ecosystem 1 has higher environmental fertility than the recipient ecosystem 2. Consumers move along-gradient, in a diffusive way that resembles how organismal movement has been represented in previous studies (e.g., Gravel et al. 2010).

Show code

# Here we set values for inorganic Nutrients input rate in each patch
# In this case, I1 > I2
I1 = 18
I2 = 2

# Leaching does not change as is still equal across ecosystems
l = 0.1

# Allocate an empty data frame to save the data in
PREDI2 <- NULL
# Allocate some rows and columns. We will replace the current numbers
# with results below Data we are calculating are (bio)mass for soil
# nutrients (N1, N2), producers (P1, P2), consumers (C1,C2),
# disperses (Q)

PREDI2 <- rbind(PREDI2, data.frame(TIME = seq(0, 100, 1), I1 = I1, I2 = I2,
    l = l, u1 = seq(0, 100, 1), u2 = seq(0, 100, 1), a1 = seq(0, 100, 1),
    a2 = seq(0, 100, 1), h1 = seq(0, 100, 1), h2 = seq(0, 100, 1), d1 = seq(0,
        100, 1), d2 = seq(0, 100, 1), g = seq(0, 100, 1), m = seq(0, 100,
        1), c = seq(0, 100, 1), e1 = seq(0, 100, 1), e2 = seq(0, 100, 1),
    N1 = seq(0, 100, 1), P1 = seq(0, 100, 1), C1 = seq(0, 100, 1), N2 = seq(0,
        100, 1), P2 = seq(0, 100, 1), C2 = seq(0, 100, 1), Q = seq(0, 100,
        1), FLUX_P1 = seq(0, 100, 1), FLUX_C1 = seq(0, 100, 1), FLUX_P2 = seq(0,
        100, 1), FLUX_C2 = seq(0, 100, 1), FLUX_Eco_1 = seq(0, 100, 1),
    FLUX_Eco_2 = seq(0, 100, 1), FLUX_Q = seq(0, 100, 1), PROD_P1 = seq(0,
        100, 1), PROD_C1 = seq(0, 100, 1), PROD_P2 = seq(0, 100, 1), PROD_C2 = seq(0,
        100, 1)))

# We want to start simulations at 0
i1 = 0

for (i in seq(0, 9999, 1)) {
    # go to the next row
    i1 = i1 + 1

    # Let's sample the parameter values from DATA. We will sample
    # these x times as defined by the i for loop above producers
    # uptake rates
    u1 = DATA1[i1, 2]  # in ecosystem 1
    u2 = DATA1[i1, 3]  # in ecosystem 2
    # consumers attack rates
    a1 = DATA1[i1, 4]
    a2 = DATA1[i1, 5]
    # death rates (= recycling rate)
    h1 = DATA1[i1, 6]  # death (=recycling) rate for producers in ecosystems 1
    h2 = DATA1[i1, 7]  # death (=recycling) rate for producers in ecosystem 2
    d1 = DATA1[i1, 8]  # death (=recycling) rate for consumers in ecosystem 1
    d2 = DATA1[i1, 9]  # death (=recycling) rate for consumers in ecosystem 2
    c = DATA1[i1, 10]  # death rate in the disperser's pool Q
    # movement rate to the Disperser's pool Q
    g = DATA1[i1, 11]
    # movement rate from the Disperser's pool Q
    m = DATA1[i1, 12]
    # efficiency of consumers
    e1 = DATA2[i1, 2]  # in ecosystem 1
    e2 = DATA2[i1, 3]  # in ecosystem 2

    PREDI2[i1, ] <- c(i1, I1, I2, l, u1, u2, a1, a2, h1, h2, d1, d2, g,
        m, c, e1, e2, ((d1 - d1 * e1 + g) * h1 + a1 * e1 * I1)/(a1 * e1 *
            l + (d1 - d1 * e1 + g) * u1), (d1 + g)/(a1 * e1), (e1 * (-h1 *
            l + I1 * u1))/(a1 * e1 * l + (d1 - d1 * e1 + g) * u1), (-a1 *
            e1 * (d2 * (-1 + e2) * h2 - a2 * e2 * I2) * l * (c + m) + d2 *
            (-1 + e2) * (d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 *
            (e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 * e2 * I2 * (c +
                m) * u1 + g * m * (h1 * l - I1 * u1))))/((c + m) * (a1 *
            e1 * l + (d1 - d1 * e1 + g) * u1) * (a2 * e2 * l - d2 * (-1 +
            e2) * u2)), (l * (-d2 * (d1 * (-1 + e1) - g) * h2 * (c + m) *
            u1 + a2 * e1 * g * m * (-h1 * l + I1 * u1)) + d2 * (d1 * e1 *
            I2 * (c + m) * u1 - (d1 + g) * I2 * (c + m) * u1 + e1 * g *
            m * (h1 * l - I1 * u1)) * u2 + a1 * d2 * e1 * l * (c + m) *
            (h2 * l - I2 * u2))/(a2 * e2 * h2 * l * (c + m) * (a1 * e1 *
            l + (d1 - d1 * e1 + g) * u1) - a2 * (a1 * e1 * e2 * I2 * l *
            (c + m) + e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 * e2 *
            I2 * (c + m) * u1 + g * m * (h1 * l - I1 * u1))) * u2), ((e1 *
            g * m * (-h1 * l + I1 * u1) * u2)/((c + m) * (a1 * e1 * l +
            (d1 - d1 * e1 + g) * u1)) + e2 * (-h2 * l + I2 * u2))/(a2 *
            e2 * l - d2 * (-1 + e2) * u2), (e1 * g * (-h1 * l + I1 * u1))/((c +
            m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1)), ((d1 + g) *
            h1)/(a1 * e1), (d1 * e1 * (-h1 * l + I1 * u1))/(a1 * e1 * l +
            (d1 - d1 * e1 + g) * u1), (h2 * (l * (-d2 * (d1 * (-1 + e1) -
            g) * h2 * (c + m) * u1 + a2 * e1 * g * m * (-h1 * l + I1 *
            u1)) + d2 * (d1 * e1 * I2 * (c + m) * u1 - (d1 + g) * I2 *
            (c + m) * u1 + e1 * g * m * (h1 * l - I1 * u1)) * u2 + a1 *
            d2 * e1 * l * (c + m) * (h2 * l - I2 * u2)))/(a2 * e2 * h2 *
            l * (c + m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1) - a2 *
            (a1 * e1 * e2 * I2 * l * (c + m) + e2 * (d1 + g) * I2 * (c +
                m) * u1 - e1 * (d1 * e2 * I2 * (c + m) * u1 + g * m * (h1 *
                l - I1 * u1))) * u2), (d2 * ((e1 * g * m * (-h1 * l + I1 *
            u1) * u2)/((c + m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1)) +
            e2 * (-h2 * l + I2 * u2)))/(a2 * e2 * l - d2 * (-1 + e2) *
            u2), ((d1 + g) * h1)/(a1 * e1) + (d1 * e1 * (-h1 * l + I1 *
            u1))/(a1 * e1 * l + (d1 - d1 * e1 + g) * u1), (h2 * (l * (-d2 *
            (d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 * e1 * g * m *
            (-h1 * l + I1 * u1)) + d2 * (d1 * e1 * I2 * (c + m) * u1 -
            (d1 + g) * I2 * (c + m) * u1 + e1 * g * m * (h1 * l - I1 *
            u1)) * u2 + a1 * d2 * e1 * l * (c + m) * (h2 * l - I2 * u2)))/(a2 *
            e2 * h2 * l * (c + m) * (a1 * e1 * l + (d1 - d1 * e1 + g) *
            u1) - a2 * (a1 * e1 * e2 * I2 * l * (c + m) + e2 * (d1 + g) *
            I2 * (c + m) * u1 - e1 * (d1 * e2 * I2 * (c + m) * u1 + g *
            m * (h1 * l - I1 * u1))) * u2) + (d2 * ((e1 * g * m * (-h1 *
            l + I1 * u1) * u2)/((c + m) * (a1 * e1 * l + (d1 - d1 * e1 +
            g) * u1)) + e2 * (-h2 * l + I2 * u2)))/(a2 * e2 * l - d2 *
            (-1 + e2) * u2), (c * e1 * g * (-h1 * l + I1 * u1))/((c + m) *
            (a1 * e1 * l + (d1 - d1 * e1 + g) * u1)), ((d1 + g) * ((d1 -
            d1 * e1 + g) * h1 + a1 * e1 * I1) * u1)/(a1 * e1 * (a1 * e1 *
            l + (d1 - d1 * e1 + g) * u1)), (e1 * (d1 + g) * (-h1 * l +
            I1 * u1))/(a1 * e1 * l + (d1 - d1 * e1 + g) * u1), ((-a1 *
            e1 * (d2 * (-1 + e2) * h2 - a2 * e2 * I2) * l * (c + m) + d2 *
            (-1 + e2) * (d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 *
            (e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 * e2 * I2 * (c +
                m) * u1 + g * m * (h1 * l - I1 * u1)))) * u2 * (l * (-d2 *
            (d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 * e1 * g * m *
            (-h1 * l + I1 * u1)) + d2 * (d1 * e1 * I2 * (c + m) * u1 -
            (d1 + g) * I2 * (c + m) * u1 + e1 * g * m * (h1 * l - I1 *
            u1)) * u2 + a1 * d2 * e1 * l * (c + m) * (h2 * l - I2 * u2)))/((c +
            m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1) * (a2 * e2 * l -
            d2 * (-1 + e2) * u2) * (a2 * e2 * h2 * l * (c + m) * (a1 *
            e1 * l + (d1 - d1 * e1 + g) * u1) - a2 * (a1 * e1 * e2 * I2 *
            l * (c + m) + e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 *
            e2 * I2 * (c + m) * u1 + g * m * (h1 * l - I1 * u1))) * u2)),
        (e2 * l * (-a1 * d2 * e1 * h2 * l * (c + m) + d2 * (d1 * (-1 +
            e1) - g) * h2 * (c + m) * u1 + a2 * e1 * g * m * (h1 * l -
            I1 * u1)) + d2 * e2 * (a1 * e1 * I2 * l * (c + m) + (d1 + g) *
            I2 * (c + m) * u1 - e1 * (d1 * I2 * (c + m) * u1 + g * m *
            (h1 * l - I1 * u1))) * u2)/((c + m) * (a1 * e1 * l + (d1 -
            d1 * e1 + g) * u1) * (a2 * e2 * l - d2 * (-1 + e2) * u2)))
}

# calculate meta-ecosystem biomass stocks for each trophic
# compartment, a.k.a., meta-ecosystem productivity
PREDI2$Ntot <- PREDI2$N1 + PREDI2$N2
PREDI2$Ptot <- PREDI2$P1 + PREDI2$P2
PREDI2$Ctot <- PREDI2$C1 + PREDI2$C2

# now calculate recycling flux for the local ecosystem
PREDI2$FLUX_Eco1_check <- PREDI2$FLUX_P1 + PREDI2$FLUX_C1
PREDI2$FLUX_Eco2_check <- PREDI2$FLUX_P2 + PREDI2$FLUX_C2

# and the meta-ecosystem recycling flux
PREDI2$FLUX_Ptot <- PREDI2$FLUX_P1 + PREDI2$FLUX_P2
PREDI2$FLUX_Ctot <- PREDI2$FLUX_C1 + PREDI2$FLUX_C2
PREDI2$MetaEcoFlux <- PREDI2$FLUX_Eco_1 + PREDI2$FLUX_Eco_2 - PREDI2$FLUX_Q

# Finally, calculate the production for each compartment at the
# meta-ecosystem scale
PREDI2 <- dplyr::mutate(PREDI2, PROD_Ptot = PROD_P1 + PROD_P2, .after = "PROD_C2")
PREDI2 <- dplyr::mutate(PREDI2, PROD_Ctot = PROD_C1 + PROD_C2, .after = "PROD_Ptot")

# Equilibrium stocks less than or equal to 0 make no biological
# sense, so we will flag them in the PREDI2 dataset.
PREDI2 <- PREDI2 %>%
    dplyr::mutate(., biosense = ifelse(N1 > 0 & P1 > 0 & C1 > 0 & N2 >
        0 & P2 > 0 & C2 > 0 & Q > 0, "yes", "no"), .after = MetaEcoFlux) %>%
    dplyr::mutate_at(vars(biosense), factor)

Stability Analyses

Here we perform stability analyses for the model in this scenario where I1 >> I2. The rationale and code are the same as above.

Show code

# create an empty dataframe
MathStabI2 <- NULL

for (i in 1:nrow(PREDI2)) {

    # fetch parameter values from the equilibrium simulations df
    I1 = PREDI2$I1[i]
    I2 = PREDI2$I2[i]
    l = PREDI2$l[i]
    u1 = PREDI2$u1[i]
    u2 = PREDI2$u2[i]
    a1 = PREDI2$a1[i]
    a2 = PREDI2$a2[i]
    h1 = PREDI2$h1[i]
    h2 = PREDI2$h2[i]
    d1 = PREDI2$d1[i]
    d2 = PREDI2$d2[i]
    g = PREDI2$g[i]
    m = PREDI2$m[i]
    c = PREDI2$c[i]
    e1 = PREDI2$e1[i]
    e2 = PREDI2$e2[i]

    # Put in the solutions to the equilibrium values here:
    dN1 = PREDI2$N1[i]
    dP1 = PREDI2$P1[i]
    dC1 = PREDI2$C1[i]
    dN2 = PREDI2$N2[i]
    dP2 = PREDI2$P2[i]
    dC2 = PREDI2$C2[i]

    # create the Jacobian from Mathematica NOTE: in transposing from
    # Mathematica, the Jacobian is here assembled by row
    Jacob <- rbind(c(-l - dP1 * u1, h1 - dN1 * u1, d1, 0, 0, 0), c(dP1 *
        u1, -a1 * dC1 - h1 + dN1 * u1, -a1 * dP1, 0, 0, 0), c(0, a1 * dC1 *
        e1, -d1 - g + a1 * e1 * dP1, 0, 0, 0), c(0, 0, 0, -l - dP2 * u2,
        h2 - dN2 * u2, d2), c(0, 0, 0, dP2 * u2, -a2 * dC2 - h2 + dN2 *
        u2, -a2 * dP2), c(0, 0, (g * m)/(c + m), 0, a2 * dC2 * e2, -d2 +
        a2 * e2 * dP2))

    # Check to see if the real part of the leading eigenvalue is less
    # than 0 or not - if it is than the system is stable.
    if (max(Re(base::eigen(Jacob)$values)) < 0) {
        stable <- "stable"
    } else {
        stable <- "unstable"
    }

    MathStabI2 <- rbind(MathStabI2, data.frame(TIME = i, I1, I2, l, u1,
        u2, a1, a2, h1, h2, d1, d2, g, m, c, e1, e2, dN1, dP1, dC1, dN2,
        dP2, dC2, EiV1 = Re(eigen(Jacob)$values[1]), EiV2 = Re(eigen(Jacob)$values[2]),
        EiV3 = Re(eigen(Jacob)$values[3]), EiV4 = Re(eigen(Jacob)$values[4]),
        EiV5 = Re(eigen(Jacob)$values[5]), EiV6 = Re(eigen(Jacob)$values[6]),
        maxEv = max(Re(base::eigen(Jacob)$values)), stable = stable, biosense = PREDI2$biosense[i]))
}

MathStabI2$stable <- as.factor(MathStabI2$stable)

# separate unstable equilibria to work with later
MathStabI2US <- subset(MathStabI2, MathStabI2$stable == "unstable")

3.24% of the stability analyses runs produced unstable results (i.e., 324 out of 10000 iterations). Now we check if these are the same ones that produce equilibria where state variables are < 0 at equilibrium (i.e., not biologicaly plausible).

Show code

I2biononsense <- subset(PREDI2, PREDI2$biosense == "no")

identical(as.numeric(MathStabI2US[, "TIME"]), I2biononsense[, "TIME"])

[1] TRUE

Since it appears that they are, let’s add the stable/unstable information to PREDI2.

Show code

PREDI2 <- left_join(PREDI2, select(MathStabI2, !c(dN1:maxEv)), by = c("TIME",
    "I1", "I2", "l", "u1", "u2", "a1", "a2", "h1", "h2", "d1", "d2", "g",
    "m", "c", "e1", "e2", "biosense"))

Now, let’s exclude the unstable parameter sets from the analyses below and from our results. Then, we will sample 1000 iterations out of PREDI2 to use in later figures.

Show code

PREDI2pos <- filter(PREDI2, PREDI2$biosense == "yes" & PREDI2$stable ==
    "stable")

# sample PREDpos to only use 1000 random simulation results
PREDI2_sample1000 <- PREDI2pos[sample(nrow(PREDI2pos), size = 1000), ]

PREDI2_1k <- droplevels(PREDI2_sample1000)

Graphs

Show code

# pivot the dataframe
PREDI2_biomass_long <- pivot_longer(PREDI2_1k, names_to = "Compartment",
    values_to = "Stock", cols = c(N1, P1, C1, N2, P2, C2, Q, Ntot, Ptot,
        Ctot)) %>%
    dplyr::mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
        "P1" | Compartment == "C1", "Donor", if_else(Compartment == "N2" |
        Compartment == "P2" | Compartment == "C2", "Recipient", if_else(Compartment ==
        "Q", "Dispersers' Pool", "Meta-ecosystem")))) %>%
    dplyr::mutate(., Ecosystem = fct_relevel(Ecosystem, "Donor", "Recipient",
        "Meta-ecosystem")) %>%
    dplyr::filter(., Ecosystem != "Dispersers' Pool")

PREDI2_biomass_long$Compartment <- as_factor(PREDI2_biomass_long$Compartment)
PREDI2_biomass_long$Ecosystem <- as_factor(PREDI2_biomass_long$Ecosystem)

comparts_medians <- PREDI2_biomass_long %>%
    dplyr::group_by(Compartment, Ecosystem) %>%
    dplyr::summarise(., median = median(Stock), .groups = "keep")

I2bm <- ggplot(data = PREDI2_biomass_long, aes(x = Compartment, y = Stock,
    col = Ecosystem)) + geom_boxjitter(alpha = 0.25, outlier.intersect = TRUE,
    outlier.shape = 25, jitter.shape = 21) + stat_summary(geom = "label",
    fun = quantile, fun.args = list(probs = 0.5), aes(label = round(after_stat(y),
        2)), position = position_dodge(width = 0.2), show.legend = F, size = 3,
    label.padding = unit(0.15, "lines"), label.size = 0.15, fontface = "bold",
    color = "black") + facet_grid(. ~ Ecosystem, scales = "free") + scale_color_manual(values = met.brewer("Egypt",
    4)) + ggtitle("Organic Biomass and Nutrients Stock") + theme_pubr() +
    coord_cartesian(y = c(0, 100)) + theme(legend.position = "none")

I2bm

Figure 6: Effects of along-gradient, diffusive consumers movement in a meta-ecosystem. Compared to the control scenario, the release of primary producers in ecosystem 1 is larger following movement of consumers out of this ecosystem (confront this figure with Figure 1). Environmental fertility values: I₁ = 18, I₂ = 2. All other specifications as in Figure 1.

Show code

FluxPREDI2 <- PREDI2_1k %>%
    dplyr::select(., TIME:c, FLUX_P1:FLUX_Eco_2, MetaEcoFlux) %>%
    tidyr::pivot_longer(., names_to = "Compartment", values_to = "Flux",
        cols = c(FLUX_P1:MetaEcoFlux)) %>%
    dplyr::mutate(., Scale = if_else(Compartment == "FLUX_P1" | Compartment ==
        "FLUX_C1" | Compartment == "FLUX_Eco_1", "Donor", if_else(Compartment ==
        "FLUX_P2" | Compartment == "FLUX_C2" | Compartment == "FLUX_Eco_2",
        "Recipient", "Meta-ecosystem"))) %>%
    dplyr::mutate(., Scale = fct_relevel(Scale, "Donor", "Recipient", "Meta-ecosystem"))

FluxPREDI2$Scale <- as.factor(FluxPREDI2$Scale)
FluxPREDI2$Compartment <- as.factor(FluxPREDI2$Compartment)

comparts_medians <- FluxPREDI2 %>%
    dplyr::group_by(Compartment, Scale) %>%
    dplyr::summarise(., median = median(Flux), .groups = "keep")

I2nf <- FluxPREDI2 %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, c("FLUX_P1",
        "FLUX_C1", "FLUX_Eco_1", "FLUX_P2", "FLUX_C2", "FLUX_Eco_2", "MetaEcoFlux"))) %>%
    ggplot(aes(x = Compartment, y = Flux)) + geom_boxjitter(outlier.intersect = TRUE,
    alpha = 0.1, outlier.shape = 25, jitter.shape = 21) + stat_summary(geom = "label",
    fun = quantile, fun.args = list(probs = 0.5), aes(label = round(after_stat(y),
        2)), position = position_dodge(width = 0.2), show.legend = F, size = 3,
    label.padding = unit(0.15, "lines"), label.size = 0.15, fontface = "bold",
    color = "black") + facet_grid(. ~ Scale, scales = "free") + coord_cartesian(y = c(0,
    300)) + scale_x_discrete(labels = c(FLUX_C1 = "C1", FLUX_P1 = "P1",
    FLUX_Eco_1 = "Donor", FLUX_C2 = "C2", FLUX_P2 = "P2", FLUX_Eco_2 = "Recipient",
    MetaEcoFlux = "Meta-ecosystem")) + ggtitle("Nutrient Flux") + theme_pubr() +
    theme(legend.position = "none", axis.title.x = element_blank())

I2nf

Figure 7: With a higher fertility in the ecosystem 1 (i.e., I₁ >> I₂) nutrient flux in this ecosystem increases, most likely as a consequence of the release primary producers. Conversely, ecosystem 2 shows lower and less variable levels of nutrient flux. At the meta-ecosystem scale, nutrient flux is higher than in either local ecosystems but lower than in the case of equal environmental fertility across ecosystems (Figure 5). Environmental fertility values: I₁ = 18, I₂ = 2. All other specifications as in Figure 5.

Show code

# I2all <- I2bm + I2nf + plot_annotation(tag_levels = 'a', tag_prefix
# = '(', tag_suffix = ')') + plot_layout(ncol = 1, nrow = 2)

# ggsave(I2all, filename = '../Results/LowI2BmNf.pdf', device =
# 'pdf', dpi = 300, width = 8, height = 6)

Against-gradient consumer movement

The code chunk below analyses the experimental scenario of higher environmental fertility in ecosystem 2. In this case, consumer movement happens counter to the resource availability gradient (i.e., non-diffusive movement; Gounand et al. (2018)).

Show code

# Here we set values for inorganic Nutrients input rate in each patch
# In this case, I1 < I2
I1 = 2
I2 = 18

# leaching rate
l = 0.1

# Allocate an empty data frame to save the data in
PREDI1 <- NULL
# Allocate some rows and columns. We will replace the current numbers
# with results below Data we are calculating are (bio)mass for soil
# nutrients (N1, N2), producers (P1, P2), consumers (C1,C2),
# disperses (Q)

PREDI1 <- rbind(PREDI1, data.frame(TIME = seq(0, 100, 1), I1 = I1, I2 = I2,
    l = l, u1 = seq(0, 100, 1), u2 = seq(0, 100, 1), a1 = seq(0, 100, 1),
    a2 = seq(0, 100, 1), h1 = seq(0, 100, 1), h2 = seq(0, 100, 1), d1 = seq(0,
        100, 1), d2 = seq(0, 100, 1), g = seq(0, 100, 1), m = seq(0, 100,
        1), c = seq(0, 100, 1), e1 = seq(0, 100, 1), e2 = seq(0, 100, 1),
    N1 = seq(0, 100, 1), P1 = seq(0, 100, 1), C1 = seq(0, 100, 1), N2 = seq(0,
        100, 1), P2 = seq(0, 100, 1), C2 = seq(0, 100, 1), Q = seq(0, 100,
        1), FLUX_P1 = seq(0, 100, 1), FLUX_C1 = seq(0, 100, 1), FLUX_P2 = seq(0,
        100, 1), FLUX_C2 = seq(0, 100, 1), FLUX_Eco_1 = seq(0, 100, 1),
    FLUX_Eco_2 = seq(0, 100, 1), FLUX_Q = seq(0, 100, 1), PROD_P1 = seq(0,
        100, 1), PROD_C1 = seq(0, 100, 1), PROD_P2 = seq(0, 100, 1), PROD_C2 = seq(0,
        100, 1)))

# We want to start simulations at 0
i1 = 0

for (i in seq(0, 9999, 1)) {
    # go to the next row
    i1 = i1 + 1

    # Let's sample the parameter values from DATA. We will sample
    # these x times as defined by the i for loop above producers
    # uptake rates
    u1 = DATA1[i1, 2]  # in ecosystem 1
    u2 = DATA1[i1, 3]  # in ecosystem 2
    # consumers attack rates
    a1 = DATA1[i1, 4]
    a2 = DATA1[i1, 5]
    # death rates (= recycling rate)
    h1 = DATA1[i1, 6]  # death (=recycling) rate for producers in ecosystems 1
    h2 = DATA1[i1, 7]  # death (=recycling) rate for producers in ecosystem 2
    d1 = DATA1[i1, 8]  # death (=recycling) rate for consumers in ecosystem 1
    d2 = DATA1[i1, 9]  # death (=recycling) rate for consumers in ecosystem 2
    c = DATA1[i1, 10]  # death rate in the disperser's pool Q
    # movement rate to the Disperser's pool Q
    g = DATA1[i1, 11]
    # movement rate from the Disperser's pool Q
    m = DATA1[i1, 12]
    # efficiency of consumers
    e1 = DATA2[i1, 2]  # in ecosystem 1
    e2 = DATA2[i1, 3]  # in ecosystem 2

    PREDI1[i1, ] <- c(i1, I1, I2, l, u1, u2, a1, a2, h1, h2, d1, d2, g,
        m, c, e1, e2, ((d1 - d1 * e1 + g) * h1 + a1 * e1 * I1)/(a1 * e1 *
            l + (d1 - d1 * e1 + g) * u1), (d1 + g)/(a1 * e1), (e1 * (-h1 *
            l + I1 * u1))/(a1 * e1 * l + (d1 - d1 * e1 + g) * u1), (-a1 *
            e1 * (d2 * (-1 + e2) * h2 - a2 * e2 * I2) * l * (c + m) + d2 *
            (-1 + e2) * (d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 *
            (e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 * e2 * I2 * (c +
                m) * u1 + g * m * (h1 * l - I1 * u1))))/((c + m) * (a1 *
            e1 * l + (d1 - d1 * e1 + g) * u1) * (a2 * e2 * l - d2 * (-1 +
            e2) * u2)), (l * (-d2 * (d1 * (-1 + e1) - g) * h2 * (c + m) *
            u1 + a2 * e1 * g * m * (-h1 * l + I1 * u1)) + d2 * (d1 * e1 *
            I2 * (c + m) * u1 - (d1 + g) * I2 * (c + m) * u1 + e1 * g *
            m * (h1 * l - I1 * u1)) * u2 + a1 * d2 * e1 * l * (c + m) *
            (h2 * l - I2 * u2))/(a2 * e2 * h2 * l * (c + m) * (a1 * e1 *
            l + (d1 - d1 * e1 + g) * u1) - a2 * (a1 * e1 * e2 * I2 * l *
            (c + m) + e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 * e2 *
            I2 * (c + m) * u1 + g * m * (h1 * l - I1 * u1))) * u2), ((e1 *
            g * m * (-h1 * l + I1 * u1) * u2)/((c + m) * (a1 * e1 * l +
            (d1 - d1 * e1 + g) * u1)) + e2 * (-h2 * l + I2 * u2))/(a2 *
            e2 * l - d2 * (-1 + e2) * u2), (e1 * g * (-h1 * l + I1 * u1))/((c +
            m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1)), ((d1 + g) *
            h1)/(a1 * e1), (d1 * e1 * (-h1 * l + I1 * u1))/(a1 * e1 * l +
            (d1 - d1 * e1 + g) * u1), (h2 * (l * (-d2 * (d1 * (-1 + e1) -
            g) * h2 * (c + m) * u1 + a2 * e1 * g * m * (-h1 * l + I1 *
            u1)) + d2 * (d1 * e1 * I2 * (c + m) * u1 - (d1 + g) * I2 *
            (c + m) * u1 + e1 * g * m * (h1 * l - I1 * u1)) * u2 + a1 *
            d2 * e1 * l * (c + m) * (h2 * l - I2 * u2)))/(a2 * e2 * h2 *
            l * (c + m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1) - a2 *
            (a1 * e1 * e2 * I2 * l * (c + m) + e2 * (d1 + g) * I2 * (c +
                m) * u1 - e1 * (d1 * e2 * I2 * (c + m) * u1 + g * m * (h1 *
                l - I1 * u1))) * u2), (d2 * ((e1 * g * m * (-h1 * l + I1 *
            u1) * u2)/((c + m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1)) +
            e2 * (-h2 * l + I2 * u2)))/(a2 * e2 * l - d2 * (-1 + e2) *
            u2), ((d1 + g) * h1)/(a1 * e1) + (d1 * e1 * (-h1 * l + I1 *
            u1))/(a1 * e1 * l + (d1 - d1 * e1 + g) * u1), (h2 * (l * (-d2 *
            (d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 * e1 * g * m *
            (-h1 * l + I1 * u1)) + d2 * (d1 * e1 * I2 * (c + m) * u1 -
            (d1 + g) * I2 * (c + m) * u1 + e1 * g * m * (h1 * l - I1 *
            u1)) * u2 + a1 * d2 * e1 * l * (c + m) * (h2 * l - I2 * u2)))/(a2 *
            e2 * h2 * l * (c + m) * (a1 * e1 * l + (d1 - d1 * e1 + g) *
            u1) - a2 * (a1 * e1 * e2 * I2 * l * (c + m) + e2 * (d1 + g) *
            I2 * (c + m) * u1 - e1 * (d1 * e2 * I2 * (c + m) * u1 + g *
            m * (h1 * l - I1 * u1))) * u2) + (d2 * ((e1 * g * m * (-h1 *
            l + I1 * u1) * u2)/((c + m) * (a1 * e1 * l + (d1 - d1 * e1 +
            g) * u1)) + e2 * (-h2 * l + I2 * u2)))/(a2 * e2 * l - d2 *
            (-1 + e2) * u2), (c * e1 * g * (-h1 * l + I1 * u1))/((c + m) *
            (a1 * e1 * l + (d1 - d1 * e1 + g) * u1)), ((d1 + g) * ((d1 -
            d1 * e1 + g) * h1 + a1 * e1 * I1) * u1)/(a1 * e1 * (a1 * e1 *
            l + (d1 - d1 * e1 + g) * u1)), (e1 * (d1 + g) * (-h1 * l +
            I1 * u1))/(a1 * e1 * l + (d1 - d1 * e1 + g) * u1), ((-a1 *
            e1 * (d2 * (-1 + e2) * h2 - a2 * e2 * I2) * l * (c + m) + d2 *
            (-1 + e2) * (d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 *
            (e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 * e2 * I2 * (c +
                m) * u1 + g * m * (h1 * l - I1 * u1)))) * u2 * (l * (-d2 *
            (d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 * e1 * g * m *
            (-h1 * l + I1 * u1)) + d2 * (d1 * e1 * I2 * (c + m) * u1 -
            (d1 + g) * I2 * (c + m) * u1 + e1 * g * m * (h1 * l - I1 *
            u1)) * u2 + a1 * d2 * e1 * l * (c + m) * (h2 * l - I2 * u2)))/((c +
            m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1) * (a2 * e2 * l -
            d2 * (-1 + e2) * u2) * (a2 * e2 * h2 * l * (c + m) * (a1 *
            e1 * l + (d1 - d1 * e1 + g) * u1) - a2 * (a1 * e1 * e2 * I2 *
            l * (c + m) + e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 *
            e2 * I2 * (c + m) * u1 + g * m * (h1 * l - I1 * u1))) * u2)),
        (e2 * l * (-a1 * d2 * e1 * h2 * l * (c + m) + d2 * (d1 * (-1 +
            e1) - g) * h2 * (c + m) * u1 + a2 * e1 * g * m * (h1 * l -
            I1 * u1)) + d2 * e2 * (a1 * e1 * I2 * l * (c + m) + (d1 + g) *
            I2 * (c + m) * u1 - e1 * (d1 * I2 * (c + m) * u1 + g * m *
            (h1 * l - I1 * u1))) * u2)/((c + m) * (a1 * e1 * l + (d1 -
            d1 * e1 + g) * u1) * (a2 * e2 * l - d2 * (-1 + e2) * u2)))
}

# calculate meta-ecosystem biomass stocks for each trophic
# compartment, a.k.a., meta-ecosystem productivity
PREDI1$Ntot <- PREDI1$N1 + PREDI1$N2
PREDI1$Ptot <- PREDI1$P1 + PREDI1$P2
PREDI1$Ctot <- PREDI1$C1 + PREDI1$C2

# now calculate recycling flux for the local ecosystem
PREDI1$FLUX_Eco1_check <- PREDI1$FLUX_P1 + PREDI1$FLUX_C1
PREDI1$FLUX_Eco2_check <- PREDI1$FLUX_P2 + PREDI1$FLUX_C2
# and the meta-ecosystem recycling flux
PREDI1$FLUX_Ptot <- PREDI1$FLUX_P1 + PREDI1$FLUX_P2
PREDI1$FLUX_Ctot <- PREDI1$FLUX_C1 + PREDI1$FLUX_C2
PREDI1$MetaEcoFlux <- PREDI1$FLUX_Eco_1 + PREDI1$FLUX_Eco_2 - PREDI1$FLUX_Q

# Finally, calculate the compartment production at meta-ecosystem
# scale
PREDI1 <- dplyr::mutate(PREDI1, PROD_Ptot = PROD_P1 + PROD_P2, .after = "PROD_C2")
PREDI1 <- dplyr::mutate(PREDI1, PROD_Ctot = PROD_C1 + PROD_C2, .after = "PROD_Ptot")

# Equilibrium stocks less than or equal to 0 make no biological
# sense, so we will remove them from the PRED dataset.

PREDI1 <- PREDI1 %>%
    dplyr::mutate(., biosense = ifelse(N1 > 0 & P1 > 0 & C1 > 0 & N2 >
        0 & P2 > 0 & C2 > 0 & Q > 0, "yes", "no"), .after = MetaEcoFlux) %>%
    dplyr::mutate_at(vars(biosense), factor)

Stability Analyses

Here we perform stability analyses for the scenario of higher environmental fertility in ecosystem 2, the recipient (i.e., I1 << I2). The analyses follow the established workflow (see above).

Show code

# create an empty dataframe
MathStabI1 <- NULL

for (i in 1:nrow(PREDI1)) {

    # fetch parameter values from the equilibrium simulations df
    I1 = PREDI1$I1[i]
    I2 = PREDI1$I2[i]
    l = PREDI1$l[i]
    u1 = PREDI1$u1[i]
    u2 = PREDI1$u2[i]
    a1 = PREDI1$a1[i]
    a2 = PREDI1$a2[i]
    h1 = PREDI1$h1[i]
    h2 = PREDI1$h2[i]
    d1 = PREDI1$d1[i]
    d2 = PREDI1$d2[i]
    g = PREDI1$g[i]
    m = PREDI1$m[i]
    c = PREDI1$c[i]
    e1 = PREDI1$e1[i]
    e2 = PREDI1$e2[i]

    # Put in the solutions to the equilibrium values here:
    dN1 = PREDI1$N1[i]
    dP1 = PREDI1$P1[i]
    dC1 = PREDI1$C1[i]
    dN2 = PREDI1$N2[i]
    dP2 = PREDI1$P2[i]
    dC2 = PREDI1$C2[i]

    # create the Jacobian from Mathematica NOTE: in transposing from
    # Mathematica, the Jacobian is here assembled by row
    Jacob <- rbind(c(-l - dP1 * u1, h1 - dN1 * u1, d1, 0, 0, 0), c(dP1 *
        u1, -a1 * dC1 - h1 + dN1 * u1, -a1 * dP1, 0, 0, 0), c(0, a1 * dC1 *
        e1, -d1 - g + a1 * e1 * dP1, 0, 0, 0), c(0, 0, 0, -l - dP2 * u2,
        h2 - dN2 * u2, d2), c(0, 0, 0, dP2 * u2, -a2 * dC2 - h2 + dN2 *
        u2, -a2 * dP2), c(0, 0, (g * m)/(c + m), 0, a2 * dC2 * e2, -d2 +
        a2 * e2 * dP2))

    # Check to see if the real part of the leading eigenvalue is less
    # than 0 or not - if it is than the system is stable.
    if (max(Re(base::eigen(Jacob)$values)) < 0) {
        stable <- "stable"
    } else {
        stable <- "unstable"
    }

    MathStabI1 <- rbind(MathStabI1, data.frame(TIME = i, I1, I2, l, u1,
        u2, a1, a2, h1, h2, d1, d2, g, m, c, e1, e2, dN1, dP1, dC1, dN2,
        dP2, dC2, EiV1 = Re(eigen(Jacob)$values[1]), EiV2 = Re(eigen(Jacob)$values[2]),
        EiV3 = Re(eigen(Jacob)$values[3]), EiV4 = Re(eigen(Jacob)$values[4]),
        EiV5 = Re(eigen(Jacob)$values[5]), EiV6 = Re(eigen(Jacob)$values[6]),
        maxEv = max(Re(base::eigen(Jacob)$values)), stable = stable, biosense = PREDI1$biosense[i]))
}

MathStabI1$stable <- as.factor(MathStabI1$stable)

# separate unstable equilibria to work with later
MathStabI1US <- subset(MathStabI1, MathStabI1$stable == "unstable")

2.85% of the stability analyses runs produced unstable results (i.e., 285 out of 10000 iterations). As always, let’s check whether these are also the parameter sets that produce state variable values at equilibrium that do not have biological meaning (i.e., <0).

Show code

I1biononsense <- subset(PREDI1, PREDI1$biosense == "no")

identical(as.numeric(MathStabI1US[, "TIME"]), I1biononsense[, "TIME"])

[1] TRUE

Now, since that proved to be the case, let’s add the stable/unstable information to PREDI1.

Show code

PREDI1 <- left_join(PREDI1, select(MathStabI1, !c(dN1:maxEv)), by = c("TIME",
    "I1", "I2", "l", "u1", "u2", "a1", "a2", "h1", "h2", "d1", "d2", "g",
    "m", "c", "e1", "e2", "biosense"))

Now, let’s exclude the unstable, biological nonsense parameter sets from the analyses below and from our results. Then, we sample 1000 random iterations from PREDI1 and store them for later use in figures and comparisons.

Show code

PREDI1pos <- filter(PREDI1, PREDI1$biosense == "yes" & PREDI1$stable ==
    "stable")

# sample PREDpos to only use 1000 random simulation results
PREDI1_sample1000 <- PREDI1pos[sample(nrow(PREDI1pos), size = 1000), ]

PREDI1_1k <- droplevels(PREDI1_sample1000)

Graphs

Show code

# pivot the
PREDI1_biomass_long <- pivot_longer(PREDI1_1k, names_to = "Compartment",
    values_to = "Stock", cols = c(N1, P1, C1, N2, P2, C2, Q, Ntot, Ptot,
        Ctot)) %>%
    dplyr::mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
        "C1" | Compartment == "P1", "Donor", if_else(Compartment == "N2" |
        Compartment == "P2" | Compartment == "C2", "Recipient", if_else(Compartment ==
        "Q", "Dispersers' Pool", "Meta-ecosystem")))) %>%
    dplyr::mutate(., Ecosystem = fct_relevel(Ecosystem, "Donor", "Recipient",
        "Meta-ecosystem")) %>%
    dplyr::filter(., Ecosystem != "Dispersers' Pool")

PREDI1_biomass_long$Compartment <- as_factor(PREDI1_biomass_long$Compartment)
PREDI1_biomass_long$Ecosystem <- as_factor(PREDI1_biomass_long$Ecosystem)

comparts_medians <- PREDI1_biomass_long %>%
    dplyr::group_by(Compartment, Ecosystem) %>%
    dplyr::summarise(., median = median(Stock), .groups = "keep")

I1bm <- ggplot(data = PREDI1_biomass_long, aes(x = Compartment, y = Stock,
    col = Ecosystem)) + geom_boxjitter(alpha = 0.25, outlier.intersect = TRUE,
    outlier.shape = 25, jitter.shape = 21) + stat_summary(geom = "label",
    fun = quantile, fun.args = list(probs = 0.5), aes(label = round(after_stat(y),
        2)), position = position_dodge(width = 0.2), show.legend = F, size = 3,
    label.padding = unit(0.15, "lines"), label.size = 0.15, fontface = "bold",
    color = "black") + facet_grid(. ~ Ecosystem, scales = "free") + scale_color_manual(values = met.brewer("Egypt",
    4)) + ggtitle("Organic Biomass and Nutrient Stock") + theme_pubr() +
    coord_cartesian(y = c(0, 100)) + theme(legend.position = "none")

I1bm

Figure 8: Movement of consumer against the environmental fertility gradient leads to a strong reduction of the nutrients stock in the donor ecosystem 1, stemming from the release of local primary producers from consumer control. In the recipient ecosystem 2, the nutrients stock increases beyond the level of enrichment seen in Figure 1, as local and immigrant consumer strongly suppress the local primary producers. All other specifications as in Figure 1. Environmental fertility values: I₁ = 2, I₂ = 18.

Following from the stronger trophic cascades, nutrient flux is also affected at both local and meta-ecosystem scales when consumer movement happens against the fertility gradient. Here, median nutrient flux in the recipient, and more fertile, ecosystem 2 is more than double the median flux in the case of diffusive consumer movement (Figure 7). Furthermore, it is also higher than in the case of equal local environmental fertility—i.e., lack of a fertility gradient (Figure 5). Likewise, median nutrient flux at the meta-ecosystem scale is higher than in either the along-gradient (Figure 7) and no-gradient (Figure 5) cases.

Show code

FluxPREDI1 <- PREDI1_1k %>%
    dplyr::select(TIME:c, FLUX_P1:FLUX_Eco_2, MetaEcoFlux) %>%
    tidyr::pivot_longer(names_to = "Compartment", values_to = "Flux", cols = c(FLUX_P1:MetaEcoFlux)) %>%
    dplyr::mutate(Scale = if_else(Compartment == "FLUX_P1" | Compartment ==
        "FLUX_C1" | Compartment == "FLUX_Eco_1", "Donor", if_else(Compartment ==
        "FLUX_P2" | Compartment == "FLUX_C2" | Compartment == "FLUX_Eco_2",
        "Recipient", "Meta-ecosystem"))) %>%
    dplyr::mutate(., Scale = fct_relevel(Scale, "Donor", "Recipient", "Meta-ecosystem"))

FluxPREDI1$Scale <- as.factor(FluxPREDI1$Scale)
FluxPREDI1$Compartment <- as.factor(FluxPREDI1$Compartment)

comparts_medians <- FluxPREDI1 %>%
    dplyr::group_by(Compartment, Scale) %>%
    dplyr::summarise(., median = median(Flux), .groups = "keep")

I1nf <- FluxPREDI1 %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, c("FLUX_P1",
        "FLUX_C1", "FLUX_Eco_1", "FLUX_P2", "FLUX_C2", "FLUX_Eco_2", "MetaEcoFlux"))) %>%
    ggplot(., aes(x = Compartment, y = Flux)) + geom_boxjitter(outlier.intersect = TRUE,
    alpha = 0.1, outlier.shape = 25, jitter.shape = 21) + stat_summary(geom = "label",
    fun = quantile, fun.args = list(probs = 0.5), aes(label = round(after_stat(y),
        2)), position = position_dodge(width = 0.2), show.legend = F, size = 3,
    label.padding = unit(0.15, "lines"), label.size = 0.15, fontface = "bold",
    color = "black") + facet_grid(. ~ Scale, scales = "free") + theme_pubr() +
    coord_cartesian(y = c(0, 400)) + scale_x_discrete(labels = c(FLUX_C1 = "C1",
    FLUX_P1 = "P1", FLUX_Eco_1 = "Donor", FLUX_C2 = "C2", FLUX_P2 = "P2",
    FLUX_Eco_2 = "Recipient", MetaEcoFlux = "Meta-ecosystem")) + ggtitle("Nutrient Flux") +
    theme_pubr() + theme(legend.position = "none", axis.title.x = element_blank())

I1nf

Figure 9: Nutrient flux more than doubles in ecosystem 2, when this is the more fertile of the two ecosystem and is also the recipient of consumer movement. In turn, this leads to an overall increase in the nutrient flux at the meta-ecosystem scale, compared to scenarios where environmental fertility is equal among local ecosystem (Figure 5) or is higher in the donor ecosystem (Figure 19). Environmental fertility values: I₁ = 2, I₂ = 18.

Show code

# I1all <- I1bm + I1nf + plot_annotation(tag_levels = 'a', tag_prefix
# = '(', tag_suffix = ')') + plot_layout(ncol = 1, nrow = 2)

# ggsave(I1all, filename = '../Results/LowI1BmNf.pdf', device =
# 'pdf', dpi = 300, width = 8, height = 6)

Summary of Stability Analyses

Here, we summarize the results of the model’s stability analyses, and collect this information to produce Table S1 shown in the Supplementary Materials. As the table shows, in all three scenario tested, only a small portion of the parameter sets used to simulate the model’s behavior at equilibrium produces unstable results. As detailed above, we exclude these unstable parameter sets from further investigations.

The table is saved in folder ../Results/ as a .tex file.

Show code

# summarise the number of stable and unstable parameter sets for each
# stability analysis run with the numerical approximation method and
# calculate the percentage of unstable parameter sets over the total
# iterations of the simulations
MStabSum <- MathStab %>%
    dplyr::group_by(., stable) %>%
    dplyr::summarise(., n = n(), .groups = "keep") %>%
    add_column(., Fertility = "Control") %>%
    pivot_wider(id_cols = Fertility, names_from = stable, values_from = n) %>%
    dplyr::mutate(., Percent = ((unstable/nrow(MathStab)) * 100))

I2MStabSum <- MathStabI2 %>%
    dplyr::group_by(., stable) %>%
    dplyr::summarise(., n = n(), .groups = "keep") %>%
    add_column(., Fertility = "Along-gradient consumer movement") %>%
    pivot_wider(id_cols = Fertility, names_from = stable, values_from = n) %>%
    dplyr::mutate(., Percent = ((unstable/nrow(MathStabI2)) * 100))

I1MStabSum <- MathStabI1 %>%
    dplyr::group_by(., stable) %>%
    dplyr::summarise(., n = n(), .groups = "keep") %>%
    add_column(., Fertility = "Against-gradient consumer movement") %>%
    pivot_wider(id_cols = Fertility, names_from = stable, values_from = n) %>%
    dplyr::mutate(., Percent = ((unstable/nrow(MathStabI1)) * 100))

# Bind the objects generated above into a single dataframe and then
# produce a table that shows the number of stable and unstable
# parameter sets for each stability analysis method

SAsummary <- bind_rows(MStabSum, I2MStabSum, I1MStabSum) %>%
    ungroup() %>%
    gt() %>%
    fmt_number(columns = 4, decimals = 2) %>%
    cols_label(., Fertility = "Fertility", stable = md("Stable"), unstable = md("Unstable"),
        Percent = md("\\%")) %>%
    tab_header(title = md("**Summary of stability analyses**, showing the number of 
             Stable and Unstable equilibria, and the percentage of Unstable 
             equilibria over the total number of iterations (n = 10000). 
             Rows are grouped by consumer movement scenario.")) %>%
    tab_style(style = list(cell_text(weight = "bold")), locations = cells_column_labels()) %>%
    tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())

# save it gtsave(SAsummary, filename =
# '../Results/StabilityAnalyses_SummaryTable.tex')

# view it
SAsummary

Summary of stability analyses, showing the number of Stable and Unstable equilibria, and the percentage of Unstable equilibria over the total number of iterations (n = 10000). Rows are grouped by consumer movement scenario.
Fertility	Stable	Unstable	%
Control	9845	155	1.55
Along-gradient consumer movement	9676	324	3.24
Against-gradient consumer movement	9715	285	2.85

Visual Deliverables—Ecosystem functions

Here we summarize the above results and produce Figures 2, 3, and 4 shown in the main text and Figures S1 to S5 in the Supplementary Materials. We refer the reader to those documents to see the actual figures; here we only present the code used to generate them.

Note that Figures 2, 3, and 4 in the main text use a log₁₀-transformed response ratio (\(LRR_X\), where \(X \in [N, P, C]\)) to assess change in nutrient stock and biomass accumulation, nutrient flux, and trophic compartment productivity in different consumer movement scenarios. Our reference for calculating the \(LRR_X\) are the results for the control scenario (Control scenario: gradient-neutral consumer movement), where environmental fertility is homogeneous between local ecosystems (i.e., I₁ = I₂ = 10).

Mathematically, the ratio’s formula is:

\(\displaystyle LRR_X = log_{10} \biggl(\frac{X_{i, I_1 \gtrless I_2}}{X_{i, I_1 = I_2}}\biggr)\)

where \(X \in [N, P, C]\) is the trophic compartment of interest, \(i\) is either the local or meta-ecosystem, and \(I_1 \gtrless I_2\) at the numerator represents the direction of fertility gradient in the experimental scenario of interest. At the denominator, \(I_1 = I_2\) represents the equal fertility conditions of the control scenario. Please, see the Numerical Analyses section of the manuscript as well as Supplementary Figure A.1 for a primer on how to interpret the LRR graphs.

All stable parameter sets

Show code

# Here we create a custom ggplot2 theme to use in the graphs below
consmov_theme <- theme_pubr() + theme(legend.position = "bottom", legend.text = element_text(size = 16),
    legend.title = element_text(size = 17), plot.title = element_text(size = 17),
    axis.title.x = element_text(size = 17, face = "bold"), axis.text.x = element_text(size = 15),
    axis.text.y = element_text(size = 15), axis.title.y = element_text(size = 17,
        face = "bold"), strip.text = element_text(size = 16))

We are working with the PRED, PREDI2, and PREDI1 datasets we produced above (sections Control scenario: gradient-neutral consumer movement and Experimental scenarios: along- and against-gradient consumer movement).

First, we will check whether there is overlap among the equilibria that are either unstable or without biological meaning across PRED, PREDI2, and PREDI1. Then, we will join them into a new object, EnvFert_res.

Show code

PREDusns <- filter(PRED, biosense == "no" | stable == "unstable")

PREDI2usns <- filter(PREDI2, biosense == "no" | stable == "unstable")

PREDI1usns <- filter(PREDI1, biosense == "no" | stable == "unstable")

summary(PREDusns[["TIME"]] %in% PREDI1usns[["TIME"]])

   Mode   FALSE    TRUE 
logical      61      94 

Show code

summary(PREDusns[["TIME"]] %in% PREDI2usns[["TIME"]])

   Mode   FALSE    TRUE 
logical      29     126 

Show code

EnvFert_res <- bind_rows(Equal = PRED, Along = PREDI2, Against = PREDI1,
    .id = "Fertility") %>%
    select(., biosense, Fertility:I2, N1:C2, Ntot:Ctot, FLUX_P1:FLUX_C2,
        FLUX_Ptot:FLUX_Ctot, PROD_P1:PROD_Ctot, u1:e2)

There is no overlap among the three datasets in terms of equilibria that are unstable and lack biological meaning. This is something that will need accounting for later on when we work with the LRR. For now, we can remove all equilibria that are either unstable or do not have biological sense as we are going to create a figure using the raw data.

Note that EnvFert_respos has the same number of rows as bmDistr created in Changes to biomass distribution in local and meta-ecosystem.

Show code

EnvFert_respos <- filter(EnvFert_res, biosense == "yes")

Untransformed results

We now transform EnvFert_respos to a long-format dataframe and add two categorical variables: Function and Scale. Function tells us which ecosystem function is measured. Scale shows instead whether the function is being measured at the local or meta-ecosystem scale.

Show code

EnvFert_long <- select(EnvFert_respos, !u1:e2) %>%
    pivot_longer(., cols = N1:PROD_Ctot, names_to = "EcoFunction", values_to = "Value") %>%
    dplyr::mutate(., Function = if_else(EcoFunction == "N1" | EcoFunction ==
        "N2" | EcoFunction == "P1" | EcoFunction == "P2" | EcoFunction ==
        "C1" | EcoFunction == "C2" | EcoFunction == "Ntot" | EcoFunction ==
        "Ptot" | EcoFunction == "Ctot", "Stock", if_else(EcoFunction ==
        "FLUX_P1" | EcoFunction == "FLUX_P2" | EcoFunction == "FLUX_C1" |
        EcoFunction == "FLUX_C2" | EcoFunction == "FLUX_Ptot" | EcoFunction ==
        "FLUX_Ctot", "Nutrient Flux", "Trophic Productivity"))) %>%
    dplyr::mutate(., Scale = if_else(EcoFunction == "N1" | EcoFunction ==
        "P1" | EcoFunction == "C1" | EcoFunction == "FLUX_P1" | EcoFunction ==
        "FLUX_C1" | EcoFunction == "PROD_P1" | EcoFunction == "PROD_C1",
        "Ecosystem 1", if_else(EcoFunction == "N2" | EcoFunction == "P2" |
            EcoFunction == "C2" | EcoFunction == "FLUX_P2" | EcoFunction ==
            "FLUX_C2" | EcoFunction == "PROD_P2" | EcoFunction == "PROD_C2",
            "Ecosystem 2", "Meta-ecosystem"))) %>%
    dplyr::mutate(., Fertility = factor(Fertility), Function = factor(Function),
        Scale = factor(Scale), EcoFunction = factor(EcoFunction))

head(EnvFert_long)

# A tibble: 6 × 9
  biosense Fertility  TIME    I1    I2 EcoFunction Value Funct…¹ Scale
  <fct>    <fct>     <dbl> <dbl> <dbl> <fct>       <dbl> <fct>   <fct>
1 yes      Equal         1    10    10 N1           3.07 Stock   Ecos…
2 yes      Equal         1    10    10 P1           1.43 Stock   Ecos…
3 yes      Equal         1    10    10 C1           1.31 Stock   Ecos…
4 yes      Equal         1    10    10 N2           3.26 Stock   Ecos…
5 yes      Equal         1    10    10 P2           1.55 Stock   Ecos…
6 yes      Equal         1    10    10 C2           2.68 Stock   Ecos…
# … with abbreviated variable name ¹​Function

Now that we have a dataset ready for plotting, we produce the graphs. The next code chunk produces Figure A.2 in the Supplementary Materials. Note that we show the code, but not the figure itself as it is presented in the Supplementary Materials.

Show code

StockSumm <- EnvFert_long %>%
    filter(., Function == "Stock") %>%
    dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, Nutrients = "N1",
        Nutrients = "N2", Nutrients = "Ntot", `Primary\n Producers` = "P1",
        `Primary\n Producers` = "P2", `Primary\n Producers` = "Ptot", Consumers = "C1",
        Consumers = "C2", Consumers = "Ctot")) %>%
    dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Nutrients",
        "Primary\n Producers", "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
        "Against")) %>%
    ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 10,
    linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
    scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Gradient-neutral",
        "Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
    scales = "free") + labs(title = "Nutrient Stock and Biomass") + ylab(" ") +
    xlab(" ") + coord_cartesian(ylim = c(0, 75))

FluxSumm <- EnvFert_long %>%
    filter(., Function == "Nutrient Flux") %>%
    dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, `Primary\n Producers` = "FLUX_P1",
        `Primary\n Producers` = "FLUX_P2", `Primary\n Producers` = "FLUX_Ptot",
        Consumers = "FLUX_C1", Consumers = "FLUX_C2", Consumers = "FLUX_Ctot")) %>%
    dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers",
        "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
        "Against")) %>%
    ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 50,
    linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
    scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Gradient-neutral",
        "Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
    scales = "free") + labs(title = "Nutrient Flux") + ylab(" ") + xlab(" ") +
    coord_cartesian(ylim = c(0, 350))

ProdSumm <- EnvFert_long %>%
    filter(., Function == "Trophic Productivity") %>%
    dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, `Primary\n Producers` = "PROD_P1",
        `Primary\n Producers` = "PROD_P2", `Primary\n Producers` = "PROD_Ptot",
        Consumers = "PROD_C1", Consumers = "PROD_C2", Consumers = "PROD_Ctot")) %>%
    dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers",
        "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
        "Against")) %>%
    ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 50,
    linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
    scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Gradient-neutral",
        "Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
    scales = "free") + labs(title = "Trophic Compartment Productivity") +
    ylab(" ") + xlab("Trophic compartment") + coord_cartesian(ylim = c(0,
    500))

FuncAll <- StockSumm + FluxSumm + ProdSumm + plot_layout(ncol = 1, nrow = 3,
    guides = "collect") + plot_annotation(tag_levels = "a", tag_prefix = "(",
    tag_suffix = ")") & consmov_theme

# FuncAll

# ggsave(FuncAll, filename = '../Results/FunctionSummary_10kN.pdf',
# device = 'pdf', dpi = 300, width = 11, height = 10)

Log-transformed Response Ratio

Now, we produce Figures 2, 3, and 4 from the main text. First, we are going to create new dataframes (PREDI2rr and PREDI1rr1) that contain only the response ratio values for nutrient stock and biomass accumulation, nutrient flux, and trophic compartment productivity of all state variables. These dataframes will also contain information about the fertility conditions under which we measure each ecosystem function. After calculating the response ratio, we will remove the unstable, biologically not meaningful equilibria from PREDI1 and PREDI2 before merging them in a new object, EnvFertRR.

In calculating the LRR, the difference in the number of equilibria that are either unstable or lack biological meaning among PRED, PREDI1, and PREDI2 matters. To account for this, let us first create two objects, I1NsUsrm and I2NsUsrm that contain the unstable, nonsensical equilibria that are found in PRED but are not included in either PREDI1 and PREDI2. We are going to remove the ones in PREDI1 and PREDI2 making use of the columns biosense and stable. Conversely, we will use I1NsUsrm and I2NsUsrm to remove instance of the response ratio being calculated using a denominator that is from an unstable, nonsensical equilibrium.

Show code

I1NsUsrm <- subset(PREDusns, !(PREDusns$TIME %in% PREDI1usns$TIME))

I2NsUsrm <- subset(PREDusns, !(PREDusns$TIME %in% PREDI2usns$TIME))

Then, let’s create two new objects, PREDI2rr and PREDI1rr, to store the response ratio values calculated by dividing ecosystem function values in PREDI2 and PREDI1 by those in PRED, respectively. We will then merge these two objects in a single one called EnvFertRR, and calculate the \(C{:}R\) biomass ratio.

Show code

# need to keep raw data in the same df as the ratio to select the >0
# cases for state variables

PREDI2rr <- PREDI2 %>%
    # we are going to batch-create new response ratio versions of
    # each variable we start with the 'stock' group of variables
dplyr::mutate(., N1rr = N1/PRED$N1, P1rr = P1/PRED$P1, C1rr = C1/PRED$C1,
    N2rr = N2/PRED$N2, P2rr = P2/PRED$P2, C2rr = C2/PRED$C2, Nrr = Ntot/PRED$Ntot,
    Prr = Ptot/PRED$Ptot, Crr = Ctot/PRED$Ctot) %>%
    # next, we work on the 'flux' group of variables
dplyr::mutate(., FLUX_P1rr = FLUX_P1/PRED$FLUX_P1, FLUX_C1rr = FLUX_C1/PRED$FLUX_C1,
    FLUX_P2rr = FLUX_P2/PRED$FLUX_P2, FLUX_C2rr = FLUX_C2/PRED$FLUX_C2,
    FLUX_Prr = FLUX_Ptot/PRED$FLUX_Ptot, FLUX_Crr = FLUX_Ctot/PRED$FLUX_Ctot) %>%
    # finally, we work on the 'productivity' group of variables
dplyr::mutate(., PROD_P1rr = PROD_P1/PRED$PROD_P1, PROD_C1rr = PROD_C1/PRED$PROD_C1,
    PROD_P2rr = PROD_P2/PRED$PROD_P2, PROD_C2rr = PROD_C2/PRED$PROD_C2,
    PROD_Prr = PROD_Ptot/PRED$PROD_Ptot, PROD_Crr = PROD_Ctot/PRED$PROD_Ctot) %>%
    # and now let's subset the dataframe to include the resp ratio
    # columns
dplyr::select(., TIME:I2, N1:Ctot, N1rr:PROD_Crr, biosense, stable, u1:e2) %>%
    dplyr::filter(., biosense == "yes" | stable == "yes") %>%
    dplyr::filter(., !(TIME %in% I2NsUsrm$TIME))

# repeat, but for PREDI1

PREDI1rr <- PREDI1 %>%
    # we are going to batch-create new response ratio versions of
    # each variable we start with the 'stock' group of variables
dplyr::mutate(., N1rr = N1/PRED$N1, P1rr = P1/PRED$P1, C1rr = C1/PRED$C1,
    N2rr = N2/PRED$N2, P2rr = P2/PRED$P2, C2rr = C2/PRED$C2, Nrr = Ntot/PRED$Ntot,
    Prr = Ptot/PRED$Ptot, Crr = Ctot/PRED$Ctot) %>%
    # next, we work on the 'flux' group of variables
dplyr::mutate(., FLUX_P1rr = FLUX_P1/PRED$FLUX_P1, FLUX_C1rr = FLUX_C1/PRED$FLUX_C1,
    FLUX_P2rr = FLUX_P2/PRED$FLUX_P2, FLUX_C2rr = FLUX_C2/PRED$FLUX_C2,
    FLUX_Prr = FLUX_Ptot/PRED$FLUX_Ptot, FLUX_Crr = FLUX_Ctot/PRED$FLUX_Ctot) %>%
    # finally, we work on the 'productivity' group of variables
dplyr::mutate(., PROD_P1rr = PROD_P1/PRED$PROD_P1, PROD_C1rr = PROD_C1/PRED$PROD_C1,
    PROD_P2rr = PROD_P2/PRED$PROD_P2, PROD_C2rr = PROD_C2/PRED$PROD_C2,
    PROD_Prr = PROD_Ptot/PRED$PROD_Ptot, PROD_Crr = PROD_Ctot/PRED$PROD_Ctot) %>%
    # and now let's subset the dataframe to include only the resp
    # ratio columns
dplyr::select(., TIME:I2, N1:Ctot, N1rr:PROD_Crr, biosense, stable, u1:e2) %>%
    dplyr::filter(., biosense == "yes" | stable == "yes") %>%
    dplyr::filter(., !(TIME %in% I1NsUsrm$TIME))

EnvFertRR <- bind_rows(Along = PREDI2rr, Against = PREDI1rr, .id = "Fertility") %>%
    dplyr::select(., Fertility, biosense, stable, TIME:I2, N1rr:PROD_Crr,
        N1:Ctot, u1:e2) %>%
    dplyr::mutate(., CR_Eco1 = C1/P1, CR_Eco2 = C2/P2, CR_MetaEco = Ctot/Ptot)
# now build the graph, first by pivoting to longer format, then
# graphing

Note that object EnvFertRR contains 19 301 and not 20 000 as would be expected from merging PREDI1 and PREDI2, because of the two rounds of removal of unstable equilibria after said merge. The first round removes those unstable, biologically not meaningful equilibria contained in PREDI1 and PREDI2, the second uses objects I1NsUsrm and I2NsUsrm to remove the ones that are contained in PRED but do not overlap with those in PREDI1 or PREDI2.

Now, we build the graphs. First, pivot the EnvFertRR dataframe from wide to long format.

Show code

EnvFertRR_long <- select(EnvFertRR, !(N1:e2)) %>%
    pivot_longer(., cols = N1rr:PROD_Crr, names_to = "Compartment", values_to = "Ratio") %>%
    dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" | Compartment ==
        "FLUX_P2rr" | Compartment == "FLUX_C1rr" | Compartment == "FLUX_C2rr" |
        Compartment == "FLUX_Prr" | Compartment == "FLUX_Crr", "Nutrient Flux",
        if_else(Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr" |
            Compartment == "PROD_P2rr" | Compartment == "PROD_C2rr" | Compartment ==
            "PROD_Prr" | Compartment == "PROD_Crr", "Trophic Productivity",
            "Stock"))) %>%
    dplyr::mutate(., Scale = if_else(Compartment == "N1rr" | Compartment ==
        "P1rr" | Compartment == "C1rr" | Compartment == "FLUX_P1rr" | Compartment ==
        "FLUX_C1rr" | Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr",
        "Ecosystem 1", if_else(Compartment == "N2rr" | Compartment == "P2rr" |
            Compartment == "C2rr" | Compartment == "FLUX_P2rr" | Compartment ==
            "FLUX_C2rr" | Compartment == "PROD_P2rr" | Compartment == "PROD_C2rr",
            "Ecosystem 2", "Meta-ecosystem"))) %>%
    dplyr::mutate(., Fertility = factor(Fertility), Function = factor(Function),
        Scale = factor(Scale), Compartment = factor(Compartment)) %>%
    dplyr::mutate(., biosense = fct_drop(biosense), stable = fct_drop(stable))

In the following code chunks, we subset and then plot each ecosystem function individually against the LRR to create Figure 2, 3, and 4 in the main text. Note, again, that here we show only the code used to generate the figures but not the figures themselves as these can be found in the main text.

The following code chunk produces Figure 2, showing changes in nutrient stocks and biomass accumulation among experimental and control scenarios.

Show code

ptHCsubset <- c("#DDAA33", "#BB5566")
ptHCsubset_fill <- c("#f8eed6", "#f1dde0")

EFstockLRR <- EnvFertRR_long %>%
    filter(., Function == "Stock") %>%
    dplyr::mutate(., Compartment = fct_recode(Compartment, Nutrients = "N1rr",
        Nutrients = "N2rr", Nutrients = "Nrr", `Primary\n Producers` = "P1rr",
        `Primary\n Producers` = "P2rr", `Primary\n Producers` = "Prr",
        Consumers = "C1rr", Consumers = "C2rr", Consumers = "Crr")) %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients",
        "Primary\n Producers", "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
    ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
    xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
    linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
    fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
    labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
    name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
    facet_grid(. ~ Scale, scales = "free") + labs(title = "Nutrient Stock and Biomass") +
    ylab("LRR") + xlab("Trophic compartment") + consmov_theme

# ggsave(EFstockLRR, filename = '../Results/Stock_LRR_10kN.pdf',
# device = 'pdf', dpi = 300, width = 10, height = 6)

The following code chunk produces Figure 3, which shows changes to nutrient flux among experimental and control scenarios.

Show code

EFfluxLRR <- EnvFertRR_long %>%
    filter(., Function == "Nutrient Flux") %>%
    dplyr::mutate(., Compartment = fct_recode(Compartment, `Primary\n Producers` = "FLUX_P1rr",
        `Primary\n Producers` = "FLUX_P2rr", `Primary\n Producers` = "FLUX_Prr",
        Consumers = "FLUX_C1rr", Consumers = "FLUX_C2rr", Consumers = "FLUX_Crr")) %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers",
        "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
    ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
    xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
    linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
    fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
    labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
    name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
    facet_grid(. ~ Scale, scales = "free") + labs(title = "Nutrient Flux") +
    ylab("LRR") + xlab("Trophic compartment") + consmov_theme

# ggsave(EFfluxLRR, filename = '../Results/Flux_LRR_10kN.pdf', device
# = 'pdf', dpi = 300, width = 10, height = 6)

And, finally, here we produce Figure 4, that shows change in the trophic compartment productivity among experimental and control scenarios.

Show code

EFprodLRR <- EnvFertRR_long %>%
    filter(., Function == "Trophic Productivity") %>%
    dplyr::mutate(., Compartment = fct_recode(Compartment, `Primary\n Producers` = "PROD_P1rr",
        `Primary\n Producers` = "PROD_P2rr", `Primary\n Producers` = "PROD_Prr",
        Consumers = "PROD_C1rr", Consumers = "PROD_C2rr", Consumers = "PROD_Crr")) %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers",
        "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
    ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
    xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
    linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
    fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
    labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
    name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
    facet_grid(. ~ Scale, scales = "free") + labs(title = "Trophic Compartment Productivity") +
    ylab("LRR") + xlab("Trophic compartment") + consmov_theme

# ggsave(EFprodLRR, filename = '../Results/Production_LRR_10kN.pdf',
# device = 'pdf', dpi = 300, width = 10, height = 6)

In the two sections below, we investigate if our results are robust to changes in biomass distribution within the system caused by moving consumers and produce Figures S2 to S5 in the Supplementary Materials.

Consumer movement in bottom-heavy systems

Here, we focus on those parameter sets that produce \(C{:}R < 1\). That is, on those parameter sets that produce classic biomass pyramids (McCauley et al. 2018).

As before, we begin with graphs for the untransformed data and then proceed to the LRR graphs.

Untransformed parameter sets

First, using the object EnvFert_CRbelow10 created above that exclude all \(C{:}R\) values above 10, we further exclude all \(C{:}R\) values above 1 and create a new object, EnvFert_CRbelow1.

Show code

EnvFert_CRbelow1 <- dplyr::mutate(EnvFert_respos, CR_Eco1 = C1/P1, CR_Eco2 = C2/P2,
    CR_MetaEco = Ctot/Ptot) %>%
    filter(., CR_Eco1 <= 1 & CR_Eco2 <= 1 & CR_MetaEco <= 1)

Now, we lengthen this new dataset for plotting with ggplot2.

Show code

EnvFert_CRbelow1_long <- select(EnvFert_CRbelow1, !u1:e2) %>%
    pivot_longer(., cols = N1:CR_MetaEco, names_to = "EcoFunction", values_to = "Value") %>%
    dplyr::mutate(., Function = if_else(EcoFunction == "N1" | EcoFunction ==
        "N2" | EcoFunction == "P1" | EcoFunction == "P2" | EcoFunction ==
        "C1" | EcoFunction == "C2" | EcoFunction == "Ntot" | EcoFunction ==
        "Ptot" | EcoFunction == "Ctot", "Stock", if_else(EcoFunction ==
        "FLUX_P1" | EcoFunction == "FLUX_P2" | EcoFunction == "FLUX_C1" |
        EcoFunction == "FLUX_C2" | EcoFunction == "FLUX_Ptot" | EcoFunction ==
        "FLUX_Ctot", "Nutrient Flux", ifelse(EcoFunction == "CR_Eco1" |
        EcoFunction == "CR_Eco2" | EcoFunction == "CR_MetaEco", "C:R",
        "Trophic Productivity")))) %>%
    dplyr::mutate(., Scale = if_else(EcoFunction == "N1" | EcoFunction ==
        "P1" | EcoFunction == "C1" | EcoFunction == "FLUX_P1" | EcoFunction ==
        "FLUX_C1" | EcoFunction == "PROD_P1" | EcoFunction == "PROD_C1" |
        EcoFunction == "CR_Eco1", "Ecosystem 1", if_else(EcoFunction ==
        "N2" | EcoFunction == "P2" | EcoFunction == "C2" | EcoFunction ==
        "FLUX_P2" | EcoFunction == "FLUX_C2" | EcoFunction == "PROD_P2" |
        EcoFunction == "PROD_C2" | EcoFunction == "CR_Eco2", "Ecosystem 2",
        "Meta-ecosystem"))) %>%
    dplyr::mutate(., Fertility = factor(Fertility), Function = factor(Function),
        Scale = factor(Scale), EcoFunction = factor(EcoFunction))

head(EnvFert_CRbelow1_long)

# A tibble: 6 × 9
  biosense Fertility  TIME    I1    I2 EcoFunct…¹  Value Funct…² Scale
  <fct>    <fct>     <dbl> <dbl> <dbl> <fct>       <dbl> <fct>   <fct>
1 yes      Equal         2    10    10 N1          0.773 Stock   Ecos…
2 yes      Equal         2    10    10 P1         16.7   Stock   Ecos…
3 yes      Equal         2    10    10 C1          0.429 Stock   Ecos…
4 yes      Equal         2    10    10 N2          0.858 Stock   Ecos…
5 yes      Equal         2    10    10 P2         12.4   Stock   Ecos…
6 yes      Equal         2    10    10 C2          3.20  Stock   Ecos…
# … with abbreviated variable names ¹​EcoFunction, ²​Function

And now we use the newly created dataframe EnvFert_CRbelow1_long to produce a figure analogous to Supplementary Figure A.2. Unlike previous figures, the following one does not appear in the Supplementary Material and is thus shown here.

Show code

StockSummCRbelow1 <- EnvFert_CRbelow1_long %>%
    filter(., Function == "Stock") %>%
    dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, Nutrients = "N1",
        Nutrients = "N2", Nutrients = "Ntot", `Primary\n Producers` = "P1",
        `Primary\n Producers` = "P2", `Primary\n Producers` = "Ptot", Consumers = "C1",
        Consumers = "C2", Consumers = "Ctot")) %>%
    dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Nutrients",
        "Primary\n Producers", "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
        "Against")) %>%
    ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 10,
    linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
    scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Equal",
        "Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
    scales = "free") + labs(title = "Nutrient Stock and Biomass") + ylab(" ") +
    xlab(" ") + coord_cartesian(ylim = c(0, 75))

FluxSummCRbelow1 <- EnvFert_CRbelow1_long %>%
    filter(., Function == "Nutrient Flux") %>%
    dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, `Primary\n Producers` = "FLUX_P1",
        `Primary\n Producers` = "FLUX_P2", `Primary\n Producers` = "FLUX_Ptot",
        Consumers = "FLUX_C1", Consumers = "FLUX_C2", Consumers = "FLUX_Ctot")) %>%
    dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers",
        "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
        "Against")) %>%
    ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 50,
    linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
    scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Equal",
        "Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
    scales = "free") + labs(title = "Nutrient Flux") + ylab(" ") + xlab(" ") +
    coord_cartesian(ylim = c(0, 350))

ProdSummCRbelow1 <- EnvFert_CRbelow1_long %>%
    filter(., Function == "Trophic Productivity") %>%
    dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, `Primary\n Producers` = "PROD_P1",
        `Primary\n Producers` = "PROD_P2", `Primary\n Producers` = "PROD_Ptot",
        Consumers = "PROD_C1", Consumers = "PROD_C2", Consumers = "PROD_Ctot")) %>%
    dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers",
        "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
        "Against")) %>%
    ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 50,
    linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
    scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Equal",
        "Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
    scales = "free") + labs(title = "Trophic Compartment Productivity") +
    ylab(" ") + xlab("Trophic Compartment") + coord_cartesian(ylim = c(0,
    500))

FuncAllCRbelow1 <- StockSummCRbelow1 + FluxSummCRbelow1 + ProdSummCRbelow1 +
    plot_layout(ncol = 1, nrow = 3, guides = "collect") + plot_annotation(tag_levels = "a",
    tag_prefix = "(", tag_suffix = ")") & consmov_theme

FuncAllCRbelow1

Figure 10: Change in nutrient stock and biomass (a), nutrient flux (b), and trophic compartment productivity (c) when consumers move in a heterogeneous meta-ecosystem, considering only parameters sets for which C:R < 1 (unstransformed data). Notably, we observe a clear quantative difference in the nutrient stock at the meta-ecosystem scale, compared with Figure A.2. However, this does not appear to cause a qualitative change in the patterns of local and meta-ecosystem functioning observe, as noted in Figure A.3. All specifications as in Figure A.2.

While the quantitative changes are more noticeable here, again in the distance between the upper and lower hinges of the boxes, there does not seem to be a marked change in the qualitative pattern observed in our results for any of the ecosystem functions of interest.

Let’s proceed to work on the LRR values.

Log-transformed Response Ratio

As we just did for the untransformed data, let’s use the previously created object EnvFertRR to select the parameter sets producing \(C{:}R < 1\) and store them in a new object, EnvFertRR_CRbelow1.

Show code

EnvFertRR_CRbelow1 <- dplyr::filter(EnvFertRR, CR_Eco1 <= 1 & CR_Eco2 <=
    1 & CR_MetaEco <= 1)

We now lengthen this newly created object.

Show code

EnvFertRR_CRbelow1_long <- select(EnvFertRR_CRbelow1, !(N1:e2)) %>%
    pivot_longer(., cols = N1rr:CR_MetaEco, names_to = "Compartment", values_to = "Ratio") %>%
    dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" | Compartment ==
        "FLUX_P2rr" | Compartment == "FLUX_C1rr" | Compartment == "FLUX_C2rr" |
        Compartment == "FLUX_Prr" | Compartment == "FLUX_Crr", "Nutrient Flux",
        if_else(Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr" |
            Compartment == "PROD_P2rr" | Compartment == "PROD_C2rr" | Compartment ==
            "PROD_Prr" | Compartment == "PROD_Crr", "Trophic Productivity",
            ifelse(Compartment == "CR_Eco1" | Compartment == "CR_Eco2" |
                Compartment == "CR_MetaEco", "C:R", "Stock")))) %>%
    dplyr::mutate(., Scale = if_else(Compartment == "N1rr" | Compartment ==
        "P1rr" | Compartment == "C1rr" | Compartment == "FLUX_P1rr" | Compartment ==
        "FLUX_C1rr" | Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr" |
        Compartment == "CR_Eco1", "Ecosystem 1", if_else(Compartment ==
        "N2rr" | Compartment == "P2rr" | Compartment == "C2rr" | Compartment ==
        "FLUX_P2rr" | Compartment == "FLUX_C2rr" | Compartment == "PROD_P2rr" |
        Compartment == "PROD_C2rr" | Compartment == "CR_Eco2", "Ecosystem 2",
        "Meta-ecosystem"))) %>%
    dplyr::mutate(., Fertility = factor(Fertility), Function = factor(Function),
        Scale = factor(Scale), Compartment = factor(Compartment)) %>%
    dplyr::mutate(., biosense = fct_drop(biosense), stable = fct_drop(stable))

And finally we produce Supplementary Figure A.3. Accordingly, we show the code here but not the figure itself.

Show code

# we make sure the y-axis of this figure has the same range of the
# 'all data' graph above by extracting the limits of the y-axis from
# that graph and using them in this one
LRRstock_ymin <- layer_scales(EFstockLRR)$y$get_limits()[1]
LRRstock_ymax <- layer_scales(EFstockLRR)$y$get_limits()[2]

EFstockLRR_CRbelow1 <- EnvFertRR_CRbelow1_long %>%
    filter(., Function == "Stock") %>%
    dplyr::mutate(., Compartment = fct_recode(Compartment, Nutrients = "N1rr",
        Nutrients = "N2rr", Nutrients = "Nrr", `Primary\n Producers` = "P1rr",
        `Primary\n Producers` = "P2rr", `Primary\n Producers` = "Prr",
        Consumers = "C1rr", Consumers = "C2rr", Consumers = "Crr")) %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients",
        "Primary\n Producers", "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
    ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
    xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
    linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
    fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
    labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
    name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
    facet_grid(. ~ Scale, scales = "free") + labs(title = "Nutrient Stock and Biomass") +
    ylab("LRR") + xlab(" ") + coord_cartesian(ylim = c(LRRstock_ymin, LRRstock_ymax))

# again, we grab the y-axis limits from the relevant all-data graph
LRRflux_ymin <- layer_scales(EFfluxLRR)$y$get_limits()[1]
LRRflux_ymax <- layer_scales(EFfluxLRR)$y$get_limits()[2]

EFfluxLRR_CRbelow1 <- EnvFertRR_CRbelow1_long %>%
    filter(., Function == "Nutrient Flux") %>%
    dplyr::mutate(., Compartment = fct_recode(Compartment, `Primary\n Producers` = "FLUX_P1rr",
        `Primary\n Producers` = "FLUX_P2rr", `Primary\n Producers` = "FLUX_Prr",
        Consumers = "FLUX_C1rr", Consumers = "FLUX_C2rr", Consumers = "FLUX_Crr")) %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers",
        "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
    ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
    xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
    linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
    fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
    labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
    name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
    facet_grid(. ~ Scale, scales = "free") + labs(title = "Nutrient Flux") +
    ylab("LRR") + xlab(" ") + coord_cartesian(ylim = c(LRRflux_ymin, LRRflux_ymax))

# again, we grab the y-axis limits from the relevant all-data graph
LRRprod_ymin <- layer_scales(EFprodLRR)$y$get_limits()[1]
LRRprod_ymax <- layer_scales(EFprodLRR)$y$get_limits()[2]

EFprodLRR_CRbelow1 <- EnvFertRR_CRbelow1_long %>%
    filter(., Function == "Trophic Productivity") %>%
    dplyr::mutate(., Compartment = fct_recode(Compartment, `Primary\n Producers` = "PROD_P1rr",
        `Primary\n Producers` = "PROD_P2rr", `Primary\n Producers` = "PROD_Prr",
        Consumers = "PROD_C1rr", Consumers = "PROD_C2rr", Consumers = "PROD_Crr")) %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers",
        "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
    ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
    xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
    linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
    fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
    labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
    name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
    facet_grid(. ~ Scale, scales = "free") + labs(title = "Trophic Compartment Productivity") +
    ylab("LRR") + xlab("Trophic compartment") + coord_cartesian(ylim = c(LRRprod_ymin,
    LRRprod_ymax))

EFallLRR_CRbelow1 <- EFstockLRR_CRbelow1 + EFfluxLRR_CRbelow1 + EFprodLRR_CRbelow1 +
    plot_layout(ncol = 1, nrow = 3, guides = "collect") + plot_annotation(tag_levels = "a",
    tag_prefix = "(", tag_suffix = ")") & consmov_theme

# EFallLRR_CRbelow1

# ggsave(EFallLRR_CRbelow1, filename =
# '../Results/FunctionSummary_LRR_CRbelow1.pdf', device = 'pdf', dpi
# = 300, width = 10, height = 10)

And here, we produce Supplementary Figure A.5, showing median parameter values for the subset of results where \(C{:}R < 1\). Again, we show the code but not the figure itself.

Show code

# table of median parameter values EnvFertRR_CRbelow1 %>% select(.,
# Fertility, u1:e2) %>% ungroup() %>% group_by(., Fertility) %>%
# pivot_longer(., cols = -Fertility, names_to = 'Parameter',
# values_to = 'value') %>% mutate(., Parameter =
# as_factor(Parameter)) %>% ungroup() %>% group_by(., Fertility,
# Parameter) %>% summarise(., across(.cols = value, .fns = list(mean
# = mean, sd = sd)), .groups = 'keep') %>% pivot_wider(., names_from
# = 'Fertility', values_from = c('value_mean', 'value_sd'))%>%
# gt(rowname_col = 'Parameter', groupname_col = 'Fertility') %>%
# fmt_number(., columns = 2:5, decimals = 2) %>% tab_spanner(label =
# md('Along'), columns = c('value_mean_Along', 'value_sd_Along')) %>%
# tab_spanner(label = md('Against'), columns =
# c('value_mean_Against', 'value_sd_Against')) %>%
# cols_label(value_mean_Against = 'mean', value_sd_Against = 'sd',
# value_mean_Along = 'mean', value_sd_Along = 'sd')

# %>% gtsave(., filename = '../Results/MeanParam_LRR_CRbelow1.rtf')

EnvFertRR_CRbelow1 %>%
    select(., Fertility, u1:e2) %>%
    ungroup() %>%
    group_by(., Fertility) %>%
    pivot_longer(., cols = -Fertility, names_to = "Parameter", values_to = "value") %>%
    mutate(., Parameter = as_factor(Parameter)) %>%
    ggplot(., aes(x = Parameter, y = value)) + geom_boxjitter(aes(fill = Parameter,
    col = Parameter), jitter.fill = NULL, alpha = 0.1) + scale_color_manual(values = met.brewer("Signac",
    13)) + scale_fill_manual(values = met.brewer("Signac", 13)) + theme_pubr(legend = "none") +
    ylab("Mean") + facet_grid(. ~ Fertility, labeller = as_labeller(c(Against = "Against-gradient Consumer movement",
    Along = "Along-gradient Consumer movement")))

Show code

# ggsave(filename = '../Results/MeanParam_LRR_CRbelow1.pdf', dpi =
# 300, device = 'pdf')

Consumer movement in top-heavy systems

Here we reproduce Figures 2, 3, and 4 from the main text and Supplementary Figure A.2 using only parameter sets that produce \(C{:}R \in (1,10)\). We exclude \(C{:}R > 10\) values, as these are likely extremely rare in real-world ecosystems. We begin, as usual, with the untransformed dataframe.

Untransformed parameter sets

First, let’s exclude from EnvFert_respos all the parameter sets producing \(C{:}R > 10\) and those producing \(C{:}R < 1\).

Show code

EnvFert_CRabove1 <- dplyr::mutate(EnvFert_respos, CR_Eco1 = C1/P1, CR_Eco2 = C2/P2,
    CR_MetaEco = Ctot/Ptot) %>%
    dplyr::filter(., CR_Eco1 <= 10 & CR_Eco2 <= 10 & CR_MetaEco <= 10) %>%
    dplyr::filter(., CR_Eco1 > 1 & CR_Eco2 > 1 & CR_MetaEco > 1)

Now we lengthen the dataset EnvFert_CRabove1 created in the code chunk above for plotting.

Show code

EnvFert_CRabove1_long <- select(EnvFert_CRabove1, !u1:e2) %>%
    pivot_longer(., cols = N1:CR_MetaEco, names_to = "EcoFunction", values_to = "Value") %>%
    dplyr::mutate(., Function = if_else(EcoFunction == "N1" | EcoFunction ==
        "N2" | EcoFunction == "P1" | EcoFunction == "P2" | EcoFunction ==
        "C1" | EcoFunction == "C2" | EcoFunction == "Ntot" | EcoFunction ==
        "Ptot" | EcoFunction == "Ctot", "Stock", if_else(EcoFunction ==
        "FLUX_P1" | EcoFunction == "FLUX_P2" | EcoFunction == "FLUX_C1" |
        EcoFunction == "FLUX_C2" | EcoFunction == "FLUX_Ptot" | EcoFunction ==
        "FLUX_Ctot", "Nutrient Flux", ifelse(EcoFunction == "CR_Eco1" |
        EcoFunction == "CR_Eco2" | EcoFunction == "CR_MetaEco", "C:R",
        "Trophic Productivity")))) %>%
    dplyr::mutate(., Scale = if_else(EcoFunction == "N1" | EcoFunction ==
        "P1" | EcoFunction == "C1" | EcoFunction == "FLUX_P1" | EcoFunction ==
        "FLUX_C1" | EcoFunction == "PROD_P1" | EcoFunction == "PROD_C1" |
        EcoFunction == "CR_Eco1", "Ecosystem 1", if_else(EcoFunction ==
        "N2" | EcoFunction == "P2" | EcoFunction == "C2" | EcoFunction ==
        "FLUX_P2" | EcoFunction == "FLUX_C2" | EcoFunction == "PROD_P2" |
        EcoFunction == "PROD_C2" | EcoFunction == "CR_Eco2", "Ecosystem 2",
        "Meta-ecosystem"))) %>%
    dplyr::mutate(., Fertility = factor(Fertility), Function = factor(Function),
        Scale = factor(Scale), EcoFunction = factor(EcoFunction))

head(EnvFert_CRabove1_long)

# A tibble: 6 × 9
  biosense Fertility  TIME    I1    I2 EcoFunct…¹  Value Funct…² Scale
  <fct>    <fct>     <dbl> <dbl> <dbl> <fct>       <dbl> <fct>   <fct>
1 yes      Equal        13    10    10 N1          7.29  Stock   Ecos…
2 yes      Equal        13    10    10 P1          0.848 Stock   Ecos…
3 yes      Equal        13    10    10 C1          3.26  Stock   Ecos…
4 yes      Equal        13    10    10 N2         15.5   Stock   Ecos…
5 yes      Equal        13    10    10 P2          0.894 Stock   Ecos…
6 yes      Equal        13    10    10 C2          7.23  Stock   Ecos…
# … with abbreviated variable names ¹​EcoFunction, ²​Function

And, finally, we reproduce Supplementary Figure A.2. Note that, again, this figure does not appear in the Supplementary Materials and is thus shown here.

Show code

StockSummCRabove1 <- EnvFert_CRabove1_long %>%
    filter(., Function == "Stock") %>%
    dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, Nutrients = "N1",
        Nutrients = "N2", Nutrients = "Ntot", `Primary\n Producers` = "P1",
        `Primary\n Producers` = "P2", `Primary\n Producers` = "Ptot", Consumers = "C1",
        Consumers = "C2", Consumers = "Ctot")) %>%
    dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Nutrients",
        "Primary\n Producers", "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
        "Against")) %>%
    ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 10,
    linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
    scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Equal",
        "Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
    scales = "free") + labs(title = "Nutrient Stock and Biomass") + ylab(" ") +
    xlab(" ") + coord_cartesian(ylim = c(0, 75))

FluxSummCRabove1 <- EnvFert_CRabove1_long %>%
    filter(., Function == "Nutrient Flux") %>%
    dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, `Primary\n Producers` = "FLUX_P1",
        `Primary\n Producers` = "FLUX_P2", `Primary\n Producers` = "FLUX_Ptot",
        Consumers = "FLUX_C1", Consumers = "FLUX_C2", Consumers = "FLUX_Ctot")) %>%
    dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers",
        "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
        "Against")) %>%
    ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 50,
    linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
    scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Equal",
        "Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
    scales = "free") + labs(title = "Nutrient Flux") + ylab(" ") + xlab(" ") +
    coord_cartesian(ylim = c(0, 350))

ProdSummCRabove1 <- EnvFert_CRabove1_long %>%
    filter(., Function == "Trophic Productivity") %>%
    dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, `Primary\n Producers` = "PROD_P1",
        `Primary\n Producers` = "PROD_P2", `Primary\n Producers` = "PROD_Ptot",
        Consumers = "PROD_C1", Consumers = "PROD_C2", Consumers = "PROD_Ctot")) %>%
    dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers",
        "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
        "Against")) %>%
    ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 50,
    linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
    scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Equal",
        "Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
    scales = "free") + labs(title = "Trophic Compartment Productivity") +
    ylab(" ") + xlab("Trophic compartment") + coord_cartesian(ylim = c(0,
    500))

FuncAllCRabove1 <- StockSummCRabove1 + FluxSummCRabove1 + ProdSummCRabove1 +
    plot_layout(ncol = 1, nrow = 3, guides = "collect") + plot_annotation(tag_levels = "a",
    tag_prefix = "(", tag_suffix = ")") & consmov_theme

FuncAllCRabove1

Figure 11: Change in nutrient stock and biomass (a), nutrient flux (b), and trophic compartment productivity (c) when consumers move in a heterogeneous meta-ecosystem, considering only parameter sets for which C:R > 1 and < 10 (unstransformed data). Here, biomass and nutrient stock results appear to more clearly point towards a trophic cascade at play at the meta-ecosystem scale, compared to Figure A.2. Patterns of nutrient flux and trophic compartment productivity appear qualitatively unchanged. All specifications as in Figure A.2.

Show code

# ggsave(FuncAllCRabove1, filename =
# '../Results/FunctionSummary_CRabove1.pdf', device = 'pdf', dpi =
# 300, width = 10, height = 10)

Compared with Supplementary Figure A.2, we see a clear change in the pattern of biomass accumulation—consistent with the inverted, top-heavy biomass pyramid expected when \(C{:}R > 1\) (McCauley et al. 2018).

Log-transformed Response Ratio

Here, we produce Supplementary Figure A.4. We use only parameter sets that produce \(C{:}R \in (1,10)\). First, we exclude all \(C{:}R\) values below 1.

Show code

EnvFertRR_CRabove1 <- dplyr::filter(EnvFertRR, CR_Eco1 <= 10 & CR_Eco2 <=
    10 & CR_MetaEco <= 10) %>%
    dplyr::filter(., CR_Eco1 > 1 & CR_Eco2 > 1 & CR_MetaEco > 1)

Then, we lengthen the dataset for plotting.

Show code

EnvFertRR_CRabove1_long <- select(EnvFertRR_CRabove1, !(N1:e2)) %>%
    pivot_longer(., cols = N1rr:CR_MetaEco, names_to = "Compartment", values_to = "Ratio") %>%
    dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" | Compartment ==
        "FLUX_P2rr" | Compartment == "FLUX_C1rr" | Compartment == "FLUX_C2rr" |
        Compartment == "FLUX_Prr" | Compartment == "FLUX_Crr", "Nutrient Flux",
        if_else(Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr" |
            Compartment == "PROD_P2rr" | Compartment == "PROD_C2rr" | Compartment ==
            "PROD_Prr" | Compartment == "PROD_Crr", "Trophic Productivity",
            ifelse(Compartment == "CR_Eco1" | Compartment == "CR_Eco2" |
                Compartment == "CR_MetaEco", "C:R", "Stock")))) %>%
    dplyr::mutate(., Scale = if_else(Compartment == "N1rr" | Compartment ==
        "P1rr" | Compartment == "C1rr" | Compartment == "FLUX_P1rr" | Compartment ==
        "FLUX_C1rr" | Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr" |
        Compartment == "CR_Eco1", "Ecosystem 1", if_else(Compartment ==
        "N2rr" | Compartment == "P2rr" | Compartment == "C2rr" | Compartment ==
        "FLUX_P2rr" | Compartment == "FLUX_C2rr" | Compartment == "PROD_P2rr" |
        Compartment == "PROD_C2rr" | Compartment == "CR_Eco2", "Ecosystem 2",
        "Meta-ecosystem"))) %>%
    dplyr::mutate(., Fertility = factor(Fertility), Function = factor(Function),
        Scale = factor(Scale), Compartment = factor(Compartment)) %>%
    dplyr::mutate(., biosense = fct_drop(biosense), stable = fct_drop(stable))

And now we produce Supplementary Figure A.4. Note that, as above, we do not show the figure itself here but only the code needed to produce it.

Show code

EFstockLRR_CRabove1 <- EnvFertRR_CRabove1_long %>%
    filter(., Function == "Stock") %>%
    dplyr::mutate(., Compartment = fct_recode(Compartment, Nutrients = "N1rr",
        Nutrients = "N2rr", Nutrients = "Nrr", `Primary\n Producers` = "P1rr",
        `Primary\n Producers` = "P2rr", `Primary\n Producers` = "Prr",
        Consumers = "C1rr", Consumers = "C2rr", Consumers = "Crr")) %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients",
        "Primary\n Producers", "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
    ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
    xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
    linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
    fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
    labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
    name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
    facet_grid(. ~ Scale, scales = "free") + labs(title = "Nutrient Stock and Biomass") +
    ylab("LRR") + xlab(" ") + coord_cartesian(ylim = c(LRRstock_ymin, LRRstock_ymax))

EFfluxLRR_CRabove1 <- EnvFertRR_CRabove1_long %>%
    filter(., Function == "Nutrient Flux") %>%
    dplyr::mutate(., Compartment = fct_recode(Compartment, `Primary\n Producers` = "FLUX_P1rr",
        `Primary\n Producers` = "FLUX_P2rr", `Primary\n Producers` = "FLUX_Prr",
        Consumers = "FLUX_C1rr", Consumers = "FLUX_C2rr", Consumers = "FLUX_Crr")) %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers",
        "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
    ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
    xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
    linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
    fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
    labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
    name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
    facet_grid(. ~ Scale, scales = "free") + labs(title = "Nutrient Flux") +
    ylab("LRR") + xlab(" ") + coord_cartesian(ylim = c(LRRflux_ymin, LRRflux_ymax))

EFprodLRR_CRabove1 <- EnvFertRR_CRabove1_long %>%
    filter(., Function == "Trophic Productivity") %>%
    dplyr::mutate(., Compartment = fct_recode(Compartment, `Primary\n Producers` = "PROD_P1rr",
        `Primary\n Producers` = "PROD_P2rr", `Primary\n Producers` = "PROD_Prr",
        Consumers = "PROD_C1rr", Consumers = "PROD_C2rr", Consumers = "PROD_Crr")) %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers",
        "Consumers")) %>%
    dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
    ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
    xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
    linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
    fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
    labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
    name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
    facet_grid(. ~ Scale, scales = "free") + labs(title = "Trophic Compartment Productivity") +
    ylab("LRR") + xlab("Trophic compartment") + coord_cartesian(ylim = c(LRRprod_ymin,
    LRRprod_ymax))

EFallLRR_CRabove1 <- EFstockLRR_CRabove1 + EFfluxLRR_CRabove1 + EFprodLRR_CRabove1 +
    plot_layout(ncol = 1, nrow = 3, guides = "collect") + plot_annotation(tag_levels = "a",
    tag_prefix = "(", tag_suffix = ")") & consmov_theme

# ggsave(EFallLRR_CRabove1, filename =
# '../Results/FunctionSummary_LRR_CRabove1.pdf', device = 'pdf', dpi
# = 300, width = 10, height = 10)

As above, we now plot the median parameter values for these parameter sets, which are shown in Supplementary Figure A.6. Accordingly, we only show the code to produce the figure and not the figure itself.

Show code

# table of median parameter sets

# EnvFertRR_CRabove1 %>% select(., Fertility, u1:e2) %>% ungroup()
# %>% group_by(., Fertility) %>% pivot_longer(., cols = -Fertility,
# names_to = 'Parameter', values_to = 'value') %>% mutate(.,
# Parameter = as_factor(Parameter)) %>% ungroup() %>% group_by(.,
# Fertility, Parameter) %>% summarise(., across(.cols = value, .fns =
# list(mean = mean, sd = sd)), .groups = 'keep') %>% pivot_wider(.,
# names_from = 'Fertility', values_from = c('value_mean',
# 'value_sd'))%>% gt(rowname_col = 'Parameter', groupname_col =
# 'Fertility') %>% fmt_number(., columns = 2:5, decimals = 2) %>%
# tab_spanner(label = md('Along'), columns = c('value_mean_Along',
# 'value_sd_Along')) %>% tab_spanner(label = md('Against'), columns =
# c('value_mean_Against', 'value_sd_Against')) %>%
# cols_label(value_mean_Against = 'mean', value_sd_Against = 'sd',
# value_mean_Along = 'mean', value_sd_Along = 'sd')
EnvFertRR_CRabove1 %>%
    select(., Fertility, u1:e2) %>%
    ungroup() %>%
    group_by(., Fertility) %>%
    pivot_longer(., cols = -Fertility, names_to = "Parameter", values_to = "value") %>%
    mutate(., Parameter = as_factor(Parameter)) %>%
    ggplot(., aes(x = Parameter, y = value)) + geom_boxjitter(aes(fill = Parameter,
    col = Parameter), jitter.fill = NULL, alpha = 0.1) + scale_color_manual(values = met.brewer("Signac",
    13)) + scale_fill_manual(values = met.brewer("Signac", 13)) + theme_pubr(legend = "none") +
    ylab("Mean") + facet_grid(. ~ Fertility, labeller = as_labeller(c(Against = "Against-gradient Consumer movement",
    Along = "Along-gradient Consumer movement")))

Show code

# ggsave(filename = '../Results/MeanParam_LRR_CRabove1.pdf', dpi =
# 300, device = 'pdf')

Summary tables—Ecosystem functions

Here, we produce summary tables for the results of the model’s iterations, using the LRR and \(C{:}R\) values we calculated to produce the summary plots above. For LRR, we report median values for each ecosystem function of interest in both local ecosystems and in the meta-ecosystem.

For ease of interpretation, median \(LRR < 0\) will be shown in shades of magenta based on how close they are to -1 and median \(LRR > 0\) will be shown in shades of green based on their closeness to 1.

First, we produce dataframes to contain the median values, one for each of stock, nutrient flux, and primary and secondary productivity.

Show code

# first, compute median of log10 ratio values

EnvFertLRR_tables <- dplyr::filter(EnvFertRR_long, Ratio >= 0) %>%
    dplyr::group_by(., Fertility, Function, Scale, Compartment) %>%
    dplyr::mutate(., LRR = log10(Ratio)) %>%
    dplyr::summarise(., min = min(LRR), median = median(LRR), max = max(LRR),
        .groups = "keep")

Now, we produce the tables, separating between along and against-gradient movement, as above.

Consumer movement changes local and meta-ecosystem functions

Biomass and Nutrient Stock

Show code

EFStockTab <- EnvFertLRR_tables %>%
    dplyr::ungroup() %>%
    dplyr::filter(., Function == "Stock") %>%
    pivot_wider(., names_from = Fertility, values_from = c(min, median,
        max)) %>%
    dplyr::mutate(., Compartment = fct_recode(Compartment, Consumers = "C1rr",
        Consumers = "C2rr", Consumers = "Crr", Producers = "P1rr", Producers = "P2rr",
        Producers = "Prr", Nutrients = "N1rr", Nutrients = "N2rr", Nutrients = "Nrr")) %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients",
        "Producers", "Consumers")) %>%
    gt(groupname_col = "Scale") %>%
    cols_hide("Function") %>%
    tab_spanner(label = md("**Along-gradient**"), columns = c("min_Along",
        "median_Along", "max_Along")) %>%
    tab_spanner(label = md("**Against-gradient**"), columns = c("min_Against",
        "median_Against", "max_Against")) %>%
    cols_label(Compartment = md("**Compartment**"), min_Along = md("*Min.*"),
        median_Along = md("*Median*"), max_Along = md("*Max.*"), min_Against = md("*Min.*"),
        median_Against = md("*Median*"), max_Against = md("*Max.*")) %>%
    fmt_number(., columns = c(4:9), decimals = 2) %>%
    data_color(., columns = c(6, 7), colors = col_numeric(palette = "PiYG",
        domain = c(-1, 1))) %>%
    tab_header(title = md("**Change in nutrient stock and organic biomass accumulation** in the meta-ecosystem when consumers move against or along an environmental fertility gradient."),
        subtitle = md("Values are expressed as median _LRR_ for 10000 iterations of the model.")) %>%
    tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())

# gtsave(EFStockTab, filename = 'EFbiostock.tex', path =
# '../Results/')

EFStockTab

Change in nutrient stock and organic biomass accumulation in the meta-ecosystem when consumers move against or along an environmental fertility gradient.
Values are expressed as median LRR for 10000 iterations of the model.
Compartment	Against-gradient			Along-gradient
	Min.	Median	Max.	Min.	Median	Max.
Ecosystem 1
Consumers	−3.16	−0.72	−0.70	0.26	0.26	2.45
Nutrients	−0.70	−0.16	0.00	0.00	0.12	0.26
Producers	0.00	0.00	0.00	0.00	0.00	0.00
Ecosystem 2
Consumers	−0.80	0.17	0.92	−1.35	−0.27	0.31
Nutrients	−0.68	0.11	0.26	−0.70	−0.14	0.25
Producers	0.00	0.04	1.46	−2.72	−0.12	0.00
Meta-ecosystem
Consumers	−1.16	0.08	0.59	−1.08	−0.09	0.57
Nutrients	−0.67	0.03	0.25	−0.68	−0.03	0.25
Producers	0.00	0.01	0.90	−0.80	−0.02	0.00

Nutrient Flux

Show code

EFFluxTab <- EnvFertLRR_tables %>%
    dplyr::ungroup() %>%
    dplyr::filter(., Function == "Nutrient Flux") %>%
    pivot_wider(., names_from = Fertility, values_from = c(min, median,
        max)) %>%
    dplyr::mutate(., Compartment = fct_recode(Compartment, Consumers = "FLUX_C1rr",
        Consumers = "FLUX_C2rr", Consumers = "FLUX_Crr", Producers = "FLUX_P1rr",
        Producers = "FLUX_P2rr", Producers = "FLUX_Prr")) %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers",
        "Consumers")) %>%
    gt(groupname_col = "Scale") %>%
    cols_hide("Function") %>%
    tab_spanner(label = md("**Along-gradient**"), columns = c("min_Along",
        "median_Along", "max_Along")) %>%
    tab_spanner(label = md("**Against-gradient**"), columns = c("min_Against",
        "median_Against", "max_Against")) %>%
    cols_label(Compartment = md("**Compartment**"), min_Along = md("*Min.*"),
        median_Along = md("*Median*"), max_Along = md("*Max.*"), min_Against = md("*Min.*"),
        median_Against = md("*Median*"), max_Against = md("*Max.*")) %>%
    fmt_number(., columns = c(4:9), decimals = 2) %>%
    data_color(., columns = c(6, 7), colors = col_numeric(palette = "PiYG",
        domain = c(-1, 1))) %>%
    tab_header(title = md("**Change in nutrient flux** in the meta-ecosystem when consumers move against or along an environmental fertility gradient."),
        subtitle = md("Values are expressed as median _LRR_ for 10000 iterations of the model.")) %>%
    tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())

# gtsave(EFFluxTab, filename = 'EFnutflux.tex', path = '../Results/')

EFFluxTab

Change in nutrient flux in the meta-ecosystem when consumers move against or along an environmental fertility gradient.
Values are expressed as median LRR for 10000 iterations of the model.
Compartment	Against-gradient			Along-gradient
	Min.	Median	Max.	Min.	Median	Max.
Ecosystem 1
Consumers	−3.16	−0.72	−0.70	0.26	0.26	2.45
Producers	0.00	0.00	0.00	0.00	0.00	0.00
Ecosystem 2
Consumers	−0.80	0.17	0.92	−1.35	−0.27	0.31
Producers	0.00	0.04	1.46	−2.72	−0.12	0.00
Meta-ecosystem
Consumers	−0.86	0.07	0.49	−1.06	−0.09	0.57
Producers	0.00	0.01	0.90	−0.86	−0.02	0.00

Trophic Productivity

Show code

EFProdTab <- EnvFertLRR_tables %>%
    dplyr::ungroup() %>%
    dplyr::filter(., Function == "Trophic Productivity") %>%
    pivot_wider(., names_from = Fertility, values_from = c(min, median,
        max)) %>%
    dplyr::mutate(., Compartment = fct_recode(Compartment, Consumers = "PROD_C1rr",
        Consumers = "PROD_C2rr", Consumers = "PROD_Crr", Producers = "PROD_P1rr",
        Producers = "PROD_P2rr", Producers = "PROD_Prr")) %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers",
        "Consumers")) %>%
    gt(groupname_col = "Scale") %>%
    cols_hide("Function") %>%
    tab_spanner(label = md("**Along-gradient**"), columns = c("min_Along",
        "median_Along", "max_Along")) %>%
    tab_spanner(label = md("**Against-gradient**"), columns = c("min_Against",
        "median_Against", "max_Against")) %>%
    cols_label(Compartment = md("**Compartment**"), min_Along = md("*Min.*"),
        median_Along = md("*Median*"), max_Along = md("*Max.*"), min_Against = md("*Min.*"),
        median_Against = md("*Median*"), max_Against = md("*Max.*")) %>%
    fmt_number(., columns = c(4:9), decimals = 2) %>%
    data_color(., columns = c(6, 7), colors = col_numeric(palette = "PiYG",
        domain = c(-1, 1))) %>%
    tab_header(title = md("**Change in trophic compartment productivity** in the meta-ecosystem when consumers move _along-_ or _against-gradient_ of nutrient availability."),
        subtitle = md("Values are expressed as median _LRR_ for 10 000 iterations of the model.")) %>%
    tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())

# gtsave(EFProdTab, filename = 'EFprod.tex', path = '../Results/')

EFProdTab

Change in trophic compartment productivity in the meta-ecosystem when consumers move along- or against-gradient of nutrient availability.
Values are expressed as median LRR for 10 000 iterations of the model.
Compartment	Against-gradient			Along-gradient
	Min.	Median	Max.	Min.	Median	Max.
Ecosystem 1
Consumers	−3.16	−0.72	−0.70	0.26	0.26	2.45
Producers	−0.70	−0.16	0.00	0.00	0.12	0.26
Ecosystem 2
Consumers	0.01	0.22	1.40	−2.63	−0.47	−0.01
Producers	0.00	0.18	1.40	−2.63	−0.34	0.00
Meta-ecosystem
Consumers	−0.88	0.03	0.47	−1.02	−0.04	0.70
Producers	−0.47	0.03	0.88	−0.81	−0.04	0.22

Changes to biomass distribution in local and meta-ecosystem

Here, we investigate whether different types of consumer movement between the two local ecosystems lead to changes in the distribution of biomass in the meta-ecosystem. Biomass distribution in local and meta-ecosystems can vary from bottom-heavy to top-heavy (McCauley et al. 2018). Bottom-heavy ecosystems present the classic biomass pyramid, with a large base of primary producers supporting an increasingly smaller biomass of primary and secondary consumers. This bottom-heavy distribution is typical of terrestrial ecosystems—albeit exceptions exist (Hatton et al. 2015; McCauley et al. 2018). Conversely, in top-heavy ecosystems, a small base of primary producers supports a large biomass of primary and secondary consumers. This inverted biomass pyramid, long seen as an exception, is increasingly recognized as widespread in certain systems—e.g., marine ecosystems (Sandin and Zgliczynski 2015; McCauley et al. 2018; Woodson, Schramski, and Joye 2018).

To investigate biomass distribution in our model, we use the Consumer to Resource biomass ratio, calculated as:

\(\displaystyle C{:}R = \frac{{C^{\ast}_i}}{{P^{\ast}_i}}\)

where the asterisk \(^{\ast}\) indicates the equilibrium biomass estimates for Consumer and Primary Producers and \(i \in [1, 2]\) indicates the ecosystem— either donor (1) or recipient (2). Values of \(C{:}R < 1\) identify ecosystems with a bottom-heavy biomass pyramid, whereas when \(C{:}R > 1\) the biomass pyramid is inverted and top-heavy (McCauley et al. 2018).

We begin by calculating the \(C{:}R\) ratio for each scenario investigated above. First, we create a new dataframe bmDistr, where we will store the required data: fertility conditions, movement rates values, efficiency rates values, state variable values at equilibrium, and ecosystem function values. We will also filter out those parameter sets that have no biological sense and are unstable as per the stability analyses above.

Show code

bmDistr <- bind_rows(Equal = PRED, Along = PREDI2, Against = PREDI1, .id = "Fertility") %>%
    select(., biosense, Fertility:e2, N1:C2, Ntot:Ctot) %>%
    filter(., biosense == "yes") %>%
    dplyr::mutate(., Eco1 = C1/P1, Eco2 = C2/P2, MetaEco = Ctot/Ptot)

Let’s look at the density distribution of these data. First, we exclude parameter sets that produce extreme values of \(C{:}R\)—that is, \(C{:}R > 10\).

Show code

ecoCR <- filter(bmDistr, Eco1 <= 10 & Eco2 <= 10 & MetaEco <= 10)

# we can save this object to speed up later work
saveRDS(ecoCR, file = "../Results/ecoCR.rds")

Now, we pivot the resulting dataset and plot the density distribution of \(C{:}R\) values.

Show code

fert.labs <- c("Against-gradient", "Along-gradient", "Gradient-neutral")
names(fert.labs) <- c("Against", "Along", "Equal")

scale.labs <- c("Ecosystem 1", "Ecosystem 2", "Meta-ecosystem")
names(scale.labs) <- c("Ecosystem 1", "Ecosystem 2", "Meta-ecosystem")

ecoCR_long <- pivot_longer(ecoCR, cols = c(Eco1, Eco2, MetaEco), values_to = "CRratio",
    names_to = "Scale") %>%
    mutate_at(., vars(Scale), factor) %>%
    mutate(., Scale = ifelse(Scale == "Eco1", "Ecosystem 1", ifelse(Scale ==
        "Eco2", "Ecosystem 2", "Meta-ecosystem")))

ggplot(ecoCR_long, aes(x = CRratio, fill = Fertility, col = Fertility)) +
    geom_density(aes(y = after_stat(density)), alpha = 0.1) + scale_color_highcontrast(reverse = T,
    name = "Consumer movement", labels = c("Against-gradient", "Along-gradient",
        "Gradient-neutral")) + scale_fill_highcontrast(reverse = T, name = "Consumer movement",
    labels = c("Against-gradient", "Along-gradient", "Gradient-neutral")) +
    facet_grid(Fertility ~ Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
        Scale = scale.labs)) + xlab("C:R ratio") + ylab("Density") + theme(legend.position = "top")

Figure 12: Density distribution when dataset is limited to C:R < 10. Note the different scales on the y-axis.

Show code

# ggsave(filename = '../Results/Density_CRratio_lessThan10.pdf',
# device = 'pdf', dpi = 300, width = 10, height = 7)

Movement from donor

Here we look at how the \(C{:}R\) values change with increasing movement from ecosystem 1, the donor, towards the dispersers’ pool Q (i.e., parameter g in the model). The graph below uses the full dataset created above, ecoCR_long, and fits linear models to each combination of Fertility scenario (gradient-neutral, along-, and against-gradient) and Scale (ecosystem 1, ecosystem 2, and meta-ecosystem).

Show code

ecoCRg_lm <- ecoCR_long %>%
    nest_by(Fertility, Scale) %>%
    mutate(model = list(lm(CRratio ~ g, data = data)))

ecoCRg_lmOutput <- ecoCRg_lm %>%
    summarise(glance(model))

ecoCRg_lmOutput$g_slope <- NA
ecoCRg_lmOutput$g_se <- NA

for (i in 1:length(ecoCRg_lm[["model"]])) {
    ecoCRg_lmOutput$g_slope[i] <- ecoCRg_lm$model[[i]]$coefficient[2]
    ecoCRg_lmOutput$g_se[i] <- summary(ecoCRg_lm$model[[i]])$coefficient["g",
        "Std. Error"]
}

ecoCRg_lmOutput %>%
    gt(groupname_col = "Scale") %>%
    fmt_number(columns = c(r.squared:p.value, logLik:deviance, g_slope:g_se),
        decimals = 3) %>%
    cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
        df.residual = "Residual df", nobs = "n", g_slope = "g Coeff.",
        g_se = "g SE") %>%
    gtsave(., filename = "../Results/CR_g_lm_res.rtf")

ecoCR_long %>%
    ggplot(., aes(x = g, y = CRratio)) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
    direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
    linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
    quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
    Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
    Scale = scale.labs)) + ylab("C:R ratio") + xlab("Movement rate from Ecosystem 1 to Dispersers' Pool (g)")

Figure 13: Distribution of C:R values lower than 10 as the movement rate from ecosystem 1 to the dispersers’ pool (g) increases. Darker shades of blue correspond to a higher count of C:R value in a given bin. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data. C:R values appear to follow a decreasing trend as g increases in Ecosystem 1 and in the Meta-ecosystem, whereas they appear stable or slightly increasing in Ecosystem 2.

Show code

# ggtitle('Binned distribution of C:R', subtitle = 'C:R values
# [0,10].')

# ggsave(filename = '../Results/CRratioVg_binned.pdf', device =
# 'pdf', dpi = 300, width = 12, height = 12)

Now, let’s explore the relationships in a bit more detail. We do this by splitting the dataframe in two portions: \(C{:}R > 1\) and \(C{:}R < 1\). Here, we take a look at the relationship for \(C{:}R\) values \(> 1\) but \(< 10\).

Show code

ecoCRg_lm_above1 <- ecoCR_long %>%
    filter(., CRratio > 1) %>%
    nest_by(Fertility, Scale) %>%
    mutate(model = list(lm(CRratio ~ g, data = data)))

ecoCRg_lmOutput_above1 <- ecoCRg_lm_above1 %>%
    summarise(glance(model))

ecoCRg_lmOutput_above1$g_slope <- NA
ecoCRg_lmOutput_above1$g_se <- NA

for (i in 1:length(ecoCRg_lm_above1[["model"]])) {
    ecoCRg_lmOutput_above1$g_slope[i] <- ecoCRg_lm_above1$model[[i]]$coefficient[2]
    ecoCRg_lmOutput_above1$g_se[i] <- summary(ecoCRg_lm_above1$model[[i]])$coefficient["g",
        "Std. Error"]
}

ecoCRg_lmOutput_above1 %>%
    gt(groupname_col = "Scale") %>%
    fmt_number(columns = c(r.squared:p.value, logLik:deviance, g_slope:g_se),
        decimals = 3) %>%
    cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
        df.residual = "Residual df", nobs = "n", g_slope = "g Coeff.",
        g_se = "g SE") %>%
    tab_header(title = md("Modeling C:R ratio vs Movement rate to Dispersers' 
                        Pool (g)"),
        subtitle = md("C:R values between (1,10].")) %>%
    gtsave(., filename = "../Results/CR_g_lm_above1.rtf")

ecoCR_long %>%
    filter(., CRratio > 1) %>%
    ggplot(., aes(x = m, y = CRratio)) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
    direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
    linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
    quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
    Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
    Scale = scale.labs)) + ylab("C:R ratio") + xlab("Movement rate from Ecosystem 1 to Dispersers' Pool (g)")

Figure 14: Distribution of C:R values comprised between (1, 10] as the movement rate from ecosystem 1 to the dispersers’ pool (g) increases. As g increases, we do not observe any trend in the values of C:R. Note that C:R > 1 are generally indicative of top-heavy, inverted biomass pyramids. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

# ggsave(filename = '../Results/CRratioVg_binned_above1.pdf', device
# = 'pdf', dpi = 300, width = 12, height = 12)

And here is the relationship for \(C{:}R\) values \(< 1\) but \(> 0\).

Show code

ecoCRg_lm_below1 <- ecoCR_long %>% 
  filter(., CRratio <= 1) %>% 
  nest_by(Fertility, Scale) %>% 
  mutate(model = list(lm(CRratio ~ g, data = data)))  

ecoCRg_lmOutput_below1 <- ecoCRg_lm_below1 %>% summarise(glance(model))

ecoCRg_lmOutput_below1$g_slope <- NA
ecoCRg_lmOutput_below1$g_se <- NA

for (i in 1:length(ecoCRg_lm_below1[["model"]])) {
  ecoCRg_lmOutput_below1$g_slope[i] <- ecoCRg_lm_below1$model[[i]]$coefficient[2]
  ecoCRg_lmOutput_below1$g_se[i] <- summary(ecoCRg_lm_above1$model[[i]])$coefficient["g", "Std. Error"]
}

ecoCRg_lmOutput_below1 %>% 
  gt(groupname_col = "Scale") %>% 
  fmt_number(., columns = c(r.squared:p.value,logLik:deviance, g_slope:g_se), 
             decimals = 3) %>% 
  cols_label(., r.squared = "R^2", 
             adj.r.squared = "Adj. R^2", 
             sigma = "SE", 
             df.residual = "Residual df", 
             nobs = "n", 
             g_slope = "g Coeff.",
             g_se = "g SE") %>% 
  tab_header(title = md("Modeling C:R ratio vs Movement rate to Dispersers' 
                        Pool (g)"), 
             subtitle = md("C:R values between [0,1].")) %>% 
  gtsave(., filename = "../Results/CR_g_lm_below1.rtf")

ecoCR_long %>% 
  filter(., CRratio <= 1) %>% 
  ggplot(., aes(x = m, y = CRratio)) + 
  geom_bin2d(bins = 70) + 
  scale_fill_met_c("Hokusai2", direction = 1) +
  geom_smooth(method = "lm", se = FALSE, col = "black",
              linetype = 2, show.legend = T, size = .5) + 
  geom_quantile(method = "rq", quantiles = .5, col = "black",
                show.legend = T) +  
 facet_grid(Fertility~Scale, scales = "free", 
            labeller = labeller(Fertility = fert.labs, 
                                Scale = scale.labs)) +  
  ylab("C:R ratio") + 
  xlab("Movement rate from Ecosystem 1 to Dispersers' Pool (g)") 

Figure 15: Distribution of C:R values [0,1] when movement rate from ecosystem 1 to the dispersers’ pool (g) increases. Again, no trend is apparent save for a seeming increase in C:R values in ecosystem 2 when consumers move along-gradient (central panel). Note, however, that this trend does not scale up to the meta-ecosystem. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

#ggsave(filename = "../Results/CRratioVg_binned_below1.pdf", device = "pdf", dpi = 300, width = 12, height = 12)

Movement to recipient

Here, we investigate how the \(C{:}R\) biomass ratio changes with increasing values of m, the movement rate of consumers leaving the dispersers’ pool Q to enter ecosystem 2. Here is the general relationship.

Show code

ecoCRm_lm <- ecoCR_long %>%
    nest_by(Fertility, Scale) %>%
    mutate(model = list(lm(CRratio ~ m, data = data)))

ecoCRm_lmOutput <- ecoCRm_lm %>%
    summarise(glance(model))

ecoCRm_lmOutput$m_slope <- NA
ecoCRm_lmOutput$m_se <- NA

for (i in 1:length(ecoCRm_lm[["model"]])) {
    ecoCRm_lmOutput$m_slope[i] <- ecoCRm_lm$model[[i]]$coefficient[2]
    ecoCRm_lmOutput$m_se[i] <- summary(ecoCRm_lm$model[[i]])$coefficient["m",
        "Std. Error"]
}

ecoCRm_lmOutput %>%
    gt(groupname_col = "Scale") %>%
    fmt_number(columns = c(r.squared:p.value, logLik:deviance, m_slope:m_se),
        decimals = 3) %>%
    cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
        df.residual = "Residual df", nobs = "n", m_slope = "m Coeff.",
        m_se = "m SE") %>%
    gtsave(., filename = "../Results/CR_m_lm_res.rtf")

ecoCR_long %>%
    ggplot(., aes(x = m, y = CRratio, )) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
    direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
    linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
    quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
    Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
    Scale = scale.labs)) + ylab("C:R ratio") + xlab("Movement rate from the Dispersers' Pool to Ecosystem 2 (m)")

Figure 16: Density distribution of C:R values [0,10] as consumer movement rate from the dispersers’ pool to ecosystem 2 (m) increases. No trend is clearly visible, aside from an apparent increase in C:R in ecosystem 2 when consumers move along-gradients (central panel). Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

# ggsave(filename = '../Results/CRratioVm_binned.pdf', device =
# 'pdf', dpi = 300, width = 12, height = 12)

Here is the relationship for \(C{:}R\) values \(> 1\) but \(< 10\).

Show code

ecoCRm_lm_above1 <- ecoCR_long %>%
    filter(., CRratio > 1) %>%
    nest_by(Fertility, Scale) %>%
    mutate(model = list(lm(CRratio ~ m, data = data)))

ecoCRm_lmOutput_above1 <- ecoCRm_lm_above1 %>%
    summarise(glance(model))

ecoCRm_lmOutput_above1$m_slope <- NA
ecoCRm_lmOutput_above1$m_se <- NA

for (i in 1:length(ecoCRm_lm_above1[["model"]])) {
    ecoCRm_lmOutput_above1$m_slope[i] <- ecoCRm_lm_above1$model[[i]]$coefficients[2]
    ecoCRm_lmOutput_above1$m_se[i] <- summary(ecoCRm_lm_above1$model[[i]])$coefficient["m",
        "Std. Error"]
}

ecoCRm_lmOutput_above1 %>%
    gt(groupname_col = "Scale") %>%
    fmt_number(columns = c(r.squared:p.value, logLik:deviance, m_slope:m_se),
        decimals = 3) %>%
    cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
        df.residual = "Residual df", nobs = "n", m_slope = "m Coeff.",
        m_se = "m SE") %>%
    tab_header(title = md("Modeling C:R ratio vs Movement rate from Dispersers' 
                        Pool (m)"),
        subtitle = md("C:R values between (1,10].")) %>%
    gtsave(., filename = "../Results/CR_m_lm_above1.rtf")

ecoCR_long %>%
    filter(., CRratio > 1) %>%
    ggplot(., aes(x = m, y = CRratio)) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
    direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
    linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
    quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
    Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
    Scale = scale.labs)) + ylab("C:R ratio") + xlab("Movement rate from the Dispersers' Pool to Ecosystem 2 (m)")

Figure 17: Distribution of C:R (1, 10] as movement rate from the dispersers’ pool to ecosystem 2 (m) increases. Focusing on only C:R values indicative of inverted, top-heavy biomass pyramids shows weakly increasing trends in the value of C:R in ecosystem 1 when consumers move against-gradient (top-left panel), and in ecosystem 2 when consumers move along-gradient (central panel). We also note an apparent decrese in C:R values in ecosystem 1 as consumers move in an homogeneous meta-ecosystem (bottom-left panel). Note, however, that none of these trends appears to influence biomass dynamics at the meta-ecosystem scale (right-most column). Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

# ggsave(filename = '../Results/CRratioVm_binned_above1.pdf', device
# = 'pdf', dpi = 300, width = 12, height = 12)

And here is the relationship for \(C{:}R\) values \(< 1\) but \(> 0\).

Show code

ecoCRm_lm_below1 <- ecoCR_long %>%
    filter(., CRratio <= 1) %>%
    nest_by(Fertility, Scale) %>%
    mutate(model = list(lm(CRratio ~ m, data = data)))

ecoCRm_lmOutput_below1 <- ecoCRm_lm_below1 %>%
    summarise(glance(model))

ecoCRm_lmOutput_below1$m_slope <- NA
ecoCRm_lmOutput_below1$m_se <- NA

for (i in 1:length(ecoCRm_lm_below1[["model"]])) {
    ecoCRm_lmOutput_below1$m_slope[i] <- ecoCRm_lm_below1$model[[i]]$coefficient[2]
    ecoCRm_lmOutput_below1$m_se[i] <- summary(ecoCRm_lm_below1$model[[i]])$coefficient["m",
        "Std. Error"]
}

ecoCRm_lmOutput_below1 %>%
    gt(groupname_col = "Scale") %>%
    fmt_number(columns = c(r.squared:p.value, logLik:deviance, m_slope:m_se),
        decimals = 3) %>%
    cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
        df.residual = "Residual df", nobs = "n", m_slope = "m Coeff.",
        m_se = "m SE") %>%
    tab_header(title = md("Modeling C:R ratio vs Movement rate from Dispersers' Pool (m)"),
        subtitle = md("C:R values between [0,1].")) %>%
    gtsave(., filename = "../Results/CR_m_lm_below1.rtf")

ecoCR_long %>%
    filter(., CRratio <= 1) %>%
    ggplot(., aes(x = m, y = CRratio)) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
    direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
    linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
    quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
    Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
    Scale = scale.labs)) + ylab("C:R ratio") + xlab("Movement rate from the Dispersers' Pool to Ecosystem 2 (m)")

Figure 18: Distribution of C:R [0,1] as consumer movement rate from the dispersers pool to ecosystem 2 (m) increases. Here, we notice apparently increasing trends in ecosystem 2 when consumers move along-gradient (central panel) and in an homogeneous meta-ecosystem (center-bottom panel). We do not see any influence of these trends at the meta-ecosystem scale. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

# ggsave(filename = '../Results/CRratioVm_binned_below1.pdf', device
# = 'pdf', dpi = 300, width = 12, height = 12)

Death rate in Dispersers’ pool

Here, we investigate how the \(C{:}R\) biomass ratio changes with increasing values of c, the death rate of consumers while in the dispersers’ pool. Consumer biomass lost while transiting through Q is lost from the meta-ecosystem and does not contribute to the recycling pathways at either local or regional scale. Hence, it is important to track how biomass distribution in the system may vary with this parameter. Here is the general relationship.

Show code

ecoCRc_lm <- ecoCR_long %>%
    nest_by(Fertility, Scale) %>%
    mutate(model = list(lm(CRratio ~ c, data = data)))

ecoCRc_lmOutput <- ecoCRc_lm %>%
    summarise(glance(model))

ecoCRc_lmOutput$c_slope <- NA
ecoCRc_lmOutput$c_se <- NA

for (i in 1:length(ecoCRc_lm[["model"]])) {
    ecoCRc_lmOutput$c_slope[i] <- ecoCRc_lm$model[[i]]$coefficient[2]
    ecoCRc_lmOutput$c_se[i] <- summary(ecoCRc_lm$model[[i]])$coefficient["c",
        "Std. Error"]
}

ecoCRc_lmOutput %>%
    gt(groupname_col = "Scale") %>%
    fmt_number(columns = c(r.squared:p.value, logLik:deviance, c_slope:c_se),
        decimals = 3) %>%
    cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
        df.residual = "Residual df", nobs = "n", c_slope = "c Coeff.",
        c_se = "c SE") %>%
    gtsave(., filename = "../Results/CR_c_lm_res.rtf")

ecoCR_long %>%
    ggplot(., aes(x = c, y = CRratio, )) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
    direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
    linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
    quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
    Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
    Scale = scale.labs)) + ylab("C:R ratio") + xlab("Consumers death rate while in the Dispersers' Pool (c)")

Figure 19: Distribution of C:R values [0,10] as the consumer death rate in the dispersers’ pool (c) increases. No trend is apparent, aside from a weak decrease in ecosystem 2 (central and center-bottom panels) that have no influence on meta-ecosystem biomass dynamics. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

# ggsave(filename = '../Results/CRratioVc_binned.pdf', device =
# 'pdf', dpi = 300, width = 12, height = 12)

Here is the relationship for \(C{:}R\) values \(> 1\) but \(< 10\).

Show code

ecoCRc_lm_above1 <- ecoCR_long %>%
    filter(., CRratio > 1) %>%
    nest_by(Fertility, Scale) %>%
    mutate(model = list(lm(CRratio ~ c, data = data)))

ecoCRc_lmOutput_above1 <- ecoCRc_lm_above1 %>%
    summarise(glance(model))

ecoCRc_lmOutput_above1$c_slope <- NA
ecoCRc_lmOutput_above1$c_se <- NA

for (i in 1:length(ecoCRc_lm_above1[["model"]])) {
    ecoCRc_lmOutput_above1$c_slope[i] <- ecoCRc_lm_above1$model[[i]]$coefficients[2]
    ecoCRc_lmOutput_above1$c_se[i] <- summary(ecoCRc_lm_above1$model[[i]])$coefficient["c",
        "Std. Error"]
}

ecoCRc_lmOutput_above1 %>%
    gt(groupname_col = "Scale") %>%
    fmt_number(columns = c(r.squared:p.value, logLik:deviance, c_slope:c_se),
        decimals = 3) %>%
    cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
        df.residual = "Residual df", nobs = "n", c_slope = "c Coeff.",
        c_se = "c SE") %>%
    tab_header(title = md("Modeling C:R ratio vs death rate while in the
                        Dispersers' Pool (c)"),
        subtitle = md("C:R values between (1,10].")) %>%
    gtsave(., filename = "../Results/CR_c_lm_above1.rtf")

ecoCR_long %>%
    filter(., CRratio > 1) %>%
    ggplot(., aes(x = c, y = CRratio)) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
    direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
    linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
    quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
    Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
    Scale = scale.labs)) + ylab("C:R ratio") + xlab("Consumer death rate while in the Dispersers' Pool (c)")

Figure 20: Distribution of C:R (1,10] as consumer death rate in the dispersers’ pool (c) increases. Regardless of the type of consumer movement, we see no trend in the distribution of C:R values that identify top-heavy, inverted biomass pyramids. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

# ggsave(filename = '../Results/CRratioVc_binned_above1.pdf', device
# = 'pdf', dpi = 300, width = 12, height = 12)

And here is the relationship for \(C{:}R\) values \(< 1\) but \(> 0\).

Show code

ecoCRc_lm_below1 <- ecoCR_long %>%
    filter(., CRratio <= 1) %>%
    nest_by(Fertility, Scale) %>%
    mutate(model = list(lm(CRratio ~ c, data = data)))

ecoCRc_lmOutput_below1 <- ecoCRc_lm_below1 %>%
    summarise(glance(model))

ecoCRc_lmOutput_below1$c_slope <- NA
ecoCRc_lmOutput_below1$c_se <- NA

for (i in 1:length(ecoCRc_lm_below1[["model"]])) {
    ecoCRc_lmOutput_below1$c_slope[i] <- ecoCRc_lm_below1$model[[i]]$coefficient[2]
    ecoCRc_lmOutput_below1$c_se[i] <- summary(ecoCRc_lm_below1$model[[i]])$coefficient["c",
        "Std. Error"]
}

ecoCRc_lmOutput_below1 %>%
    gt(groupname_col = "Scale") %>%
    fmt_number(columns = c(r.squared:p.value, logLik:deviance, c_slope:c_se),
        decimals = 3) %>%
    cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
        df.residual = "Residual df", nobs = "n", c_slope = "c Coeff.",
        c_se = "c SE") %>%
    tab_header(title = md("Modeling C:R ratio vs death rate while in the 
                        Dispersers' pool (c)"),
        subtitle = md("C:R values between [0,1].")) %>%
    gtsave(., filename = "../Results/CR_c_lm_below1.rtf")

ecoCR_long %>%
    filter(., CRratio <= 1) %>%
    ggplot(., aes(x = c, y = CRratio)) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
    direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
    linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
    quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
    Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
    Scale = scale.labs)) + ylab("C:R ratio") + xlab("Consumer death rate while in the Dispersers' Pool (c)")

Figure 21: Distribution of C:R [0,1] as consumer death rate in the dispersers’ pool (c) increases. Focusing on C:R < 1 brings out a weak, negative trend in ecosystem 2 when consumers move along-gradient (central panel) and in an homogeneous meta-ecosystem (center-bottom panel). Neither of these trends have any influence at the meta-ecosystem scale. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

# ggsave(filename = '../Results/CRratioVc_binned_below1.pdf', device
# = 'pdf', dpi = 300, width = 12, height = 12)

Visual Deliverable—Biomass pyramids

Here we produce graphs of the biomass pyramids produced by consumer movement across ecosystems. These graphs allow for easy visualization of the changes in the distribution of biomass in either local ecosystem, or in the meta-ecosystem as a whole, following different types of consumer movement.

We produce four graphs, one each for the following four cases: (i) classic, bottom-heavy pyramid in the donor ecosystem 1; (ii) classic, bottom-heavy pyramid in the recipient ecosystem 2; (iii) inverted, top-heavy pyramid in the donor ecosystem 1; (iv) inverted, top-heavy pyramid in the recipient ecosystem 2. We separate between bottom- and top-heavy pyramids in either the donor or recipient ecosystems as these local conditions can influence the effects of consumer movement on biomass distribution and thus ecosystem functioning. Note that these four biomass-informed categories will return when we perform a global sensitivity analyses on our results below (Sensitivity Analyses).

The code we use here was initially developed by Dr. A. Z. Andis Arietta.

Classic donor

Show code

# first, produce a df that contains the values that will be printed
# on top of each bar in the biomass pyramids
summary_df <- ecoCR %>%
    ungroup() %>%
    select(biosense:I2, N1:Ctot, Eco1) %>%
    filter(., Eco1 < 1) %>%
    pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
    mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
        "P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
        "N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
        "Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
        "N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
        ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
            "Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
        "Nutrients", "Primary\nProducers", "Consumers")) %>%
    group_by(Fertility, Ecosystem, Compartment)

# then, produce the plot
ecoCR %>%
    ungroup() %>%
    select(biosense:I2, N1:Ctot, Eco1) %>%
    filter(Eco1 < 1) %>%
    pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
    mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
        "P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
        "N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
        "Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
        "N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
        ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
            "Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
        "Nutrients", "Primary\nProducers", "Consumers")) %>%
    ggplot(aes(x = Compartment, y = Biomass)) + geom_bar(aes(fill = Compartment,
    color = Compartment), stat = "identity", width = 0.9) + geom_bar(data = . %>%
    mutate(Biomass = -Biomass), aes(fill = Compartment, color = Compartment),
    stat = "identity", width = 0.9) + stat_summary(data = summary_df, geom = "label",
    fun.args = list(mult = 1, na.rm = T), aes(label = round(after_stat(y),
        2)), position = position_dodge(width = 1), show.legend = F, size = 2,
    label.padding = unit(0.1, "lines"), label.size = 0.1, color = "black") +
    scale_color_met_d(name = "Isfahan2", direction = 1) + scale_fill_met_d(name = "Isfahan2",
    direction = 1) + facet_grid(Fertility ~ Ecosystem, scales = "free") +
    coord_flip() + theme_classic() + theme(text = element_text(size = 14),
    axis.title = element_blank(), axis.ticks = element_blank(), axis.text = element_blank(),
    legend.position = "top")

Figure 22: Biomass pyramids for ecosystem 1 (donor), ecosystem 2 (recipient), and the meta-ecosystem when the Consumer:Resource ratio in ecosystem 1 is < 1. Bottom-heavy biomass pyramids dominate the system in all movement scenarios. However, note how against-gradient consumer movement consistently allows for—on average—higher levels of Primary Producers biomass compared to both along-gradient and gradient-neutral consumer movement. The labels inside the bars show the average biomass for each ecosystem compartment.

Show code

ggsave(filename = "../Results/Bpyram_CRbelow1_eco1.pdf", dpi = 300, device = "pdf")

Classic recipient

Show code

# first, produce a df that contains the values that will be printed
# on top of each bar in the biomass pyramids
summary_df <- ecoCR %>%
    ungroup() %>%
    select(biosense:I2, N1:Ctot, Eco2) %>%
    filter(Eco2 < 1) %>%
    pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
    mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
        "P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
        "N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
        "Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
        "N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
        ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
            "Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
        "Nutrients", "Primary\nProducers", "Consumers")) %>%
    group_by(Fertility, Ecosystem, Compartment)

# then, produce the plot
ecoCR %>%
    ungroup() %>%
    select(biosense:I2, N1:Ctot, Eco2) %>%
    filter(Eco2 < 1) %>%
    pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
    mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
        "P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
        "N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
        "Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
        "N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
        ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
            "Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
        "Nutrients", "Primary\nProducers", "Consumers")) %>%
    ggplot(aes(x = Compartment, y = Biomass)) + geom_bar(aes(fill = Compartment,
    color = Compartment), stat = "identity", width = 0.9) + geom_bar(data = . %>%
    mutate(Biomass = -Biomass), aes(fill = Compartment, color = Compartment),
    stat = "identity", width = 0.9) + stat_summary(data = summary_df, geom = "label",
    fun.args = list(mult = 1, na.rm = T), aes(label = round(after_stat(y),
        2)), position = position_dodge(width = 1), show.legend = F, size = 2,
    label.padding = unit(0.1, "lines"), label.size = 0.1, color = "black") +
    scale_color_met_d(name = "Isfahan2", direction = 1) + scale_fill_met_d(name = "Isfahan2",
    direction = 1) + facet_grid(Fertility ~ Ecosystem, scales = "free") +
    coord_flip() + theme_classic() + theme(text = element_text(size = 14),
    axis.title = element_blank(), axis.ticks = element_blank(), axis.text = element_blank(),
    legend.position = "top")

Figure 23: Biomass pyramids for ecosystem 1 (donor), ecosystem 2 (recipient), and the meta-ecosystem when the Consumer:Resource ratio in ecosystem 2 is < 1. In this case, the biomass distribution in the system is more nuanced and variable, and depends on the direction of consumer movement. Against-gradient consumer movement leads to comparable average biomass of Primary Producers in the two ecosystems. Conversely, in the along-gradient and gradient-neutral scenarios, the distribution of Primary Producers’ biomass is skewed towards ecosystem 1 (the donor). Note also how the distribution of Nutrients varies among donor and recipient ecosystems across movement scenarios.

Show code

ggsave(filename = "../Results/Bpyram_CRbelow1_eco2.pdf", dpi = 300, device = "pdf")

Inverted donor

Show code

# first, produce a df that contains the values that will be printed
# on top of each bar in the biomass pyramids
summary_df <- ecoCR %>%
    ungroup() %>%
    select(biosense:I2, N1:Ctot, Eco1) %>%
    filter(Eco1 > 1) %>%
    pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
    mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
        "P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
        "N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
        "Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
        "N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
        ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
            "Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
        "Nutrients", "Primary\nProducers", "Consumers")) %>%
    group_by(Fertility, Ecosystem, Compartment)

# then, produce the plot
ecoCR %>%
    ungroup() %>%
    select(biosense:I2, N1:Ctot, Eco1) %>%
    filter(Eco1 > 1) %>%
    pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
    mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
        "P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
        "N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
        "Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
        "N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
        ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
            "Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
        "Nutrients", "Primary\nProducers", "Consumers")) %>%
    ggplot(aes(x = Compartment, y = Biomass)) + geom_bar(aes(fill = Compartment,
    color = Compartment), stat = "identity", width = 0.9) + geom_bar(data = . %>%
    mutate(Biomass = -Biomass), aes(fill = Compartment, color = Compartment),
    stat = "identity", width = 0.9) + stat_summary(data = summary_df, geom = "label",
    fun.args = list(mult = 1, na.rm = T), aes(label = round(after_stat(y),
        2)), position = position_dodge(width = 1), show.legend = F, size = 2,
    label.padding = unit(0.1, "lines"), label.size = 0.1, color = "black") +
    scale_color_met_d(name = "Isfahan2", direction = 1) + scale_fill_met_d(name = "Isfahan2",
    direction = 1) + facet_grid(Fertility ~ Ecosystem, scales = "free") +
    coord_flip() + theme_classic() + theme(text = element_text(size = 14),
    axis.title = element_blank(), axis.ticks = element_blank(), axis.text = element_blank(),
    legend.position = "top")

Figure 24: Biomass pyramids for ecosystem 1 (donor), ecosystem 2 (recipient), and the meta-ecosystem when the Consumer:Resource ratio in ecosystem 1 is > 1. In this case, while we see top-heavy, inverted biomass pyramids in the donor ecosystem, bottom-heavy pyramids still dominate the recipient ecosystem and ultimately the meta-ecosystem as a whole. In the against-gradient scenario, we can see an accumulation of nutrients in the recipient ecosystem (resource pooling effect), whereas the along-gradient scenario leads to a more equal distribution of nutrients in the system that resembles what we find in the gradient-neutral scenario.

Show code

ggsave(filename = "../Results/Bpyram_CRabove1_eco1.pdf", dpi = 300, device = "pdf")

Inverted recipient

Show code

# first, produce a df that contains the values that will be printed
# on top of each bar in the biomass pyramids
summary_df <- ecoCR %>%
    ungroup() %>%
    select(biosense:I2, N1:Ctot, Eco2) %>%
    filter(Eco2 > 1) %>%
    pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
    mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
        "P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
        "N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
        "Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
        "N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
        ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
            "Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
        "Nutrients", "Primary\nProducers", "Consumers")) %>%
    group_by(Fertility, Ecosystem, Compartment)

# then, produce the plot
ecoCR %>%
    ungroup() %>%
    select(biosense:I2, N1:Ctot, Eco2) %>%
    filter(Eco2 > 1) %>%
    pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
    mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
        "P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
        "N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
        "Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
        "N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
        ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
            "Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
        "Nutrients", "Primary\nProducers", "Consumers")) %>%
    ggplot(aes(x = Compartment, y = Biomass)) + geom_bar(aes(fill = Compartment,
    color = Compartment), stat = "identity", width = 0.9) + geom_bar(data = . %>%
    mutate(Biomass = -Biomass), aes(fill = Compartment, color = Compartment),
    stat = "identity", width = 0.9) + stat_summary(data = summary_df, geom = "label",
    fun.args = list(mult = 1, na.rm = T), aes(label = round(after_stat(y),
        2)), position = position_dodge(width = 1), show.legend = F, size = 2,
    label.padding = unit(0.1, "lines"), label.size = 0.1, color = "black") +
    scale_color_met_d(name = "Isfahan2", direction = 1) + scale_fill_met_d(name = "Isfahan2",
    direction = 1) + facet_grid(Fertility ~ Ecosystem, scales = "free") +
    coord_flip() + theme_classic() + theme(text = element_text(size = 14),
    axis.title = element_blank(), axis.ticks = element_blank(), axis.text = element_blank(),
    legend.position = "top")

Figure 25: Biomass pyramids for ecosystem 1 (donor), ecosystem 2 (recipient), and the meta-ecosystem when the Consumer:Resource ratio in ecosystem 2 is > 1. Here, we observe the spatial trophic cascade in action as a lack of top-down control on Primary Producers leads to a highly skewed accumulation of biomass for this compartment in ecosystem 1, the donor, whereas ecosystem 2 shows much reduced autothroph presence. Conversely, both Consumer biomass and Nutrient stock are higher in ecosystem 2, the recipient, in all movement scenarios. Notice how, when the recipient ecosystem shows a top-heavy biomass pyramid, along-gradient consumer movement is the only scenario in which the distribution of nutrients is somewhat homogeneous across the landscape.

Show code

ggsave(filename = "../Results/Bpyram_CRabove1_eco2.pdf", dpi = 300, device = "pdf")

Summary tables—Biomass distribution

The three tabs below contain tables showing summary statistics for the linear models fitted to the distribution of biomass in the meta-ecosystem. The table shown in tab [Movement from ecosystem 1] is intended as a companion to Figures 13, 14, and 15. Similarly, the table in tab [Movement towards ecosystem 2] is a companion to Figures 16, 17, and 18. Finally, the table in tab Death rate in Dispersers’ pool provides details to accompany Figures 19, 20, and 21.

Movement from donor

Show code

# first, compute median of log10 ratio values

ecoCR_g_table <- bind_rows(`above 1` = ecoCRg_lmOutput_above1, `below 1` = ecoCRg_lmOutput_below1,
    .id = "CRslice") %>%
    ungroup() %>%
    select(., CRslice:adj.r.squared, g_slope:g_se) %>%
    pivot_wider(., names_from = "CRslice", values_from = c(r.squared, adj.r.squared,
        g_slope, g_se)) %>%
    gt(., groupname_col = "Scale") %>%
    tab_spanner(label = md("**C:R (1, 10]**"), columns = c("r.squared_above 1",
        "adj.r.squared_above 1", "g_slope_above 1", "g_se_above 1")) %>%
    tab_spanner(label = md("**C:R [0, 1]**"), columns = c("r.squared_below 1",
        "adj.r.squared_below 1", "g_slope_below 1", "g_se_below 1")) %>%
    fmt_number(., columns = 3:10, decimals = 3) %>%
    cols_label(., `r.squared_above 1` = md("R sq."), `adj.r.squared_above 1` = md("Adj. R sq."),
        `g_slope_above 1` = md("g Coeff."), `g_se_above 1` = md("g SE"),
        `r.squared_below 1` = md("R sq."), `adj.r.squared_below 1` = md("Adj. R sq."),
        `g_slope_below 1` = md("g Coeff."), `g_se_below 1` = md("g SE")) %>%
    # data_color(., columns = c(3:10), colors = col_numeric(palette =
    # 'PiYG', domain = c(-1,1))) %>%
tab_header(title = md("**Relationships between biomass _C:R_ and movement from ecosystem 1.**"),
    subtitle = md("The table shows the slope value (g Coeff.) and standard error (g SE) for linear regression models fitted to the biomass _C:R_ data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between _C:R_  values larger than 1 but below 10, and _C:R_ values below 1 but above 0.")) %>%
    tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())

ecoCR_g_table

Relationships between biomass C:R and movement from ecosystem 1.
The table shows the slope value (g Coeff.) and standard error (g SE) for linear regression models fitted to the biomass C:R data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between C:R values larger than 1 but below 10, and C:R values below 1 but above 0.
Fertility	C:R (1, 10]				C:R [0, 1]
	R sq.	Adj. R sq.	g Coeff.	g SE	R sq.	Adj. R sq.	g Coeff.	g SE
Ecosystem 1
Against	0.083	0.078	−0.728	0.180	0.119	0.119	−0.017	0.180
Along	0.148	0.147	−0.343	0.022	0.010	0.010	−0.009	0.022
Equal	0.177	0.176	−0.537	0.040	0.039	0.039	−0.016	0.040
Ecosystem 2
Against	0.001	0.001	−0.025	0.016	0.000	0.000	0.001	0.016
Along	0.000	0.000	0.018	0.015	0.006	0.005	0.007	0.015
Equal	0.001	0.001	−0.027	0.015	0.001	0.000	0.003	0.015
Meta-ecosystem
Against	0.034	0.033	−0.078	0.012	0.002	0.001	−0.004	0.012
Along	0.048	0.047	−0.100	0.011	0.002	0.002	−0.004	0.011
Equal	0.049	0.048	−0.106	0.012	0.005	0.005	−0.006	0.012

Movement to recipient

Show code

# first, compute median of log10 ratio values

ecoCR_m_table <- bind_rows(`above 1` = ecoCRm_lmOutput_above1, `below 1` = ecoCRm_lmOutput_below1,
    .id = "CRslice") %>%
    ungroup() %>%
    select(., CRslice:adj.r.squared, m_slope:m_se) %>%
    pivot_wider(., names_from = "CRslice", values_from = c(r.squared, adj.r.squared,
        m_slope, m_se)) %>%
    gt(., groupname_col = "Scale") %>%
    tab_spanner(label = md("**C:R (1, 10]**"), columns = c("r.squared_above 1",
        "adj.r.squared_above 1", "m_slope_above 1", "m_se_above 1")) %>%
    tab_spanner(label = md("**C:R [0, 1]**"), columns = c("r.squared_below 1",
        "adj.r.squared_below 1", "m_slope_below 1", "m_se_below 1")) %>%
    fmt_number(., columns = 3:10, decimals = 3) %>%
    cols_label(., `r.squared_above 1` = md("R sq."), `adj.r.squared_above 1` = md("Adj. R sq."),
        `m_slope_above 1` = md("m Coeff."), `m_se_above 1` = md("m SE"),
        `r.squared_below 1` = md("R sq."), `adj.r.squared_below 1` = md("Adj. R sq."),
        `m_slope_below 1` = md("m Coeff."), `m_se_below 1` = md("m SE")) %>%
    # data_color(., columns = c(3:10), colors = col_numeric(palette =
    # 'PiYG', domain = c(-,))) %>%
tab_header(title = md("**Relationships between biomass _C:R_ and movement into ecosystem 2.**"),
    subtitle = md("The table shows the slope value (m Coeff.) and standard error (m SE) for linear regression models fitted to the biomass _C:R_ data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between _C:R_  values larger than 1 but below 10, and _C:R_ values below 1 but above 0.")) %>%
    tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())

ecoCR_m_table

Relationships between biomass C:R and movement into ecosystem 2.
The table shows the slope value (m Coeff.) and standard error (m SE) for linear regression models fitted to the biomass C:R data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between C:R values larger than 1 but below 10, and C:R values below 1 but above 0.
Fertility	C:R (1, 10]				C:R [0, 1]
	R sq.	Adj. R sq.	m Coeff.	m SE	R sq.	Adj. R sq.	m Coeff.	m SE
Ecosystem 1
Against	0.003	−0.002	0.036	0.048	0.000	0.000	0.000	0.001
Along	0.000	0.000	−0.014	0.018	0.001	0.001	−0.002	0.001
Equal	0.002	0.001	−0.031	0.026	0.000	0.000	−0.001	0.001
Ecosystem 2
Against	0.000	0.000	−0.008	0.016	0.000	0.000	0.000	0.002
Along	0.001	0.001	0.026	0.015	0.016	0.016	0.012	0.001
Equal	0.000	0.000	−0.001	0.015	0.006	0.006	0.007	0.001
Meta-ecosystem
Against	0.000	−0.001	0.000	0.013	0.000	0.000	0.000	0.001
Along	0.000	0.000	0.010	0.011	0.000	0.000	0.001	0.001
Equal	0.001	0.000	−0.014	0.013	0.000	0.000	0.001	0.001

Death rate in Dispersers’ pool

Show code

# first, compute median of log10 ratio values

ecoCR_c_table <- bind_rows(`above 1` = ecoCRc_lmOutput_above1, `below 1` = ecoCRc_lmOutput_below1,
    .id = "CRslice") %>%
    ungroup() %>%
    select(., CRslice:adj.r.squared, c_slope:c_se) %>%
    pivot_wider(., names_from = "CRslice", values_from = c(r.squared, adj.r.squared,
        c_slope, c_se)) %>%
    gt(., groupname_col = "Scale") %>%
    tab_spanner(label = md("**C:R (1, 10]**"), columns = c("r.squared_above 1",
        "adj.r.squared_above 1", "c_slope_above 1", "c_se_above 1")) %>%
    tab_spanner(label = md("**C:R [0, 1]**"), columns = c("r.squared_below 1",
        "adj.r.squared_below 1", "c_slope_below 1", "c_se_below 1")) %>%
    fmt_number(., columns = 3:10, decimals = 3) %>%
    cols_label(., `r.squared_above 1` = md("R sq."), `adj.r.squared_above 1` = md("Adj. R sq."),
        `c_slope_above 1` = md("c Coeff."), `c_se_above 1` = md("c SE"),
        `r.squared_below 1` = md("R sq."), `adj.r.squared_below 1` = md("Adj. R sq."),
        `c_slope_below 1` = md("c Coeff."), `c_se_below 1` = md("c SE")) %>%
    # data_color(., columns = c(3:10), colors = col_numeric(palette =
    # 'PiYG', domain = c(-,))) %>%
tab_header(title = md("**Relationships between biomass _C:R_ and death rate in dispersers' pool.**"),
    subtitle = md("The table shows the slope value (c Coeff.) and standard error (c SE) for linear regression models fitted to the biomass _C:R_ data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between _C:R_  values larger than 1 but below 10, and _C:R_ values below 1 but above 0.")) %>%
    tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())

ecoCR_c_table

Relationships between biomass C:R and death rate in dispersers' pool.
The table shows the slope value (c Coeff.) and standard error (c SE) for linear regression models fitted to the biomass C:R data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between C:R values larger than 1 but below 10, and C:R values below 1 but above 0.
Fertility	C:R (1, 10]				C:R [0, 1]
	R sq.	Adj. R sq.	c Coeff.	c SE	R sq.	Adj. R sq.	c Coeff.	c SE
Ecosystem 1
Against	0.005	−0.001	0.047	0.050	0.000	0.000	−0.001	0.001
Along	0.000	0.000	−0.012	0.019	0.000	0.000	0.002	0.001
Equal	0.000	−0.001	0.000	0.026	0.000	0.000	0.001	0.001
Ecosystem 2
Against	0.000	0.000	0.010	0.016	0.000	0.000	0.000	0.002
Along	0.000	0.000	−0.004	0.014	0.002	0.002	−0.004	0.001
Equal	0.001	0.001	−0.030	0.015	0.002	0.002	−0.005	0.002
Meta-ecosystem
Against	0.003	0.003	0.025	0.012	0.000	0.000	0.000	0.001
Along	0.000	−0.001	−0.003	0.011	0.001	0.001	−0.003	0.001
Equal	0.000	−0.001	0.004	0.012	0.001	0.001	−0.003	0.001

Summary tables—Comparison across C:R values

Here, we produce summary tables that show how the effects of consumer movement on local and meta-ecosystem functions vary as the biomass distribution changes in the system as a result of the consumer movement itself.

We begin by collating together the required dataframes.

Show code

CRcomparison <- bind_rows(all = EnvFertRR, above_1 = EnvFertRR_CRabove1,
    below_1 = EnvFertRR_CRbelow1, .id = "CR") %>%
    select(., !(N1:Ctot)) %>%
    filter(., biosense == "yes" | stable == "stable") %>%
    mutate(., CR = as_factor(CR)) %>%
    droplevels()

Then, we lengthen the dataset as this is the preferred format for producing tables with package gt.

Show code

# the first step in the pipe remove the by-ecosystem CR data, as
# these are now store in the CR column
CRcomparison_long <- select(CRcomparison, !c(biosense, stable, CR_Eco1,
    CR_Eco2, CR_MetaEco)) %>%
    pivot_longer(., cols = N1rr:PROD_Crr, names_to = "Compartment", values_to = "Ratio") %>%
    dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" | Compartment ==
        "FLUX_P2rr" | Compartment == "FLUX_C1rr" | Compartment == "FLUX_C2rr" |
        Compartment == "FLUX_Prr" | Compartment == "FLUX_Crr", "Nutrient Flux",
        if_else(Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr" |
            Compartment == "PROD_P2rr" | Compartment == "PROD_C2rr" | Compartment ==
            "PROD_Prr" | Compartment == "PROD_Crr", "Trophic Productivity",
            "Stock"))) %>%
    dplyr::mutate(., Scale = if_else(Compartment == "N1rr" | Compartment ==
        "P1rr" | Compartment == "C1rr" | Compartment == "FLUX_P1rr" | Compartment ==
        "FLUX_C1rr" | Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr",
        "Ecosystem 1", if_else(Compartment == "N2rr" | Compartment == "P2rr" |
            Compartment == "C2rr" | Compartment == "FLUX_P2rr" | Compartment ==
            "FLUX_C2rr" | Compartment == "PROD_P2rr" | Compartment == "PROD_C2rr",
            "Ecosystem 2", "Meta-ecosystem"))) %>%
    dplyr::mutate(., Fertility = factor(Fertility), Function = factor(Function),
        Scale = factor(Scale), Compartment = factor(Compartment))

We now remove all values of the Response Ratio that are \(< 0\) and \(log_{10}\)-transform them before calculating the median.

Show code

CRcomparison_tables <- dplyr::filter(CRcomparison_long, Ratio >= 0) %>%
    dplyr::group_by(., CR, Fertility, Function, Scale, Compartment) %>%
    dplyr::mutate(., LRR = log10(Ratio)) %>%
    dplyr::summarise(., median = median(LRR), .groups = "keep")

# can add min(LRR) and max(LRR) is needed

Finally, we produce summary tables that show the differences in median LRR when considering all stable parameter sets, parameter sets that produce \(C{:}R > 1\), and parameter sets that produce \(C{:}R < 1\). As before, we produce one table per ecosystem function.

Biomass and Nutrient Accumulation

Show code

CRsummStock <- CRcomparison_tables %>%
    ungroup() %>%
    filter(., Function == "Stock") %>%
    dplyr::mutate(., Compartment = fct_recode(Compartment, Consumers = "C1rr",
        Consumers = "C2rr", Consumers = "Crr", Producers = "P1rr", Producers = "P2rr",
        Producers = "Prr", Nutrients = "N1rr", Nutrients = "N2rr", Nutrients = "Nrr")) %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients",
        "Producers", "Consumers")) %>%
    pivot_wider(., names_from = c(CR, Fertility), values_from = median) %>%
    gt(., groupname_col = "Scale") %>%
    fmt_number(., columns = 4:9, decimals = 2) %>%
    tab_spanner(label = md("**Along-gradient**"), columns = c("all_Along",
        "above_1_Along", "below_1_Along")) %>%
    tab_spanner(label = md("**Against-gradient**"), columns = c("all_Against",
        "above_1_Against", "below_1_Against")) %>%
    cols_label(Compartment = md("**Compartment**"), all_Against = md("*All C:R*"),
        above_1_Against = md("*C:R > 1*"), below_1_Against = md("*C:R < 1*"),
        all_Along = md("*All C:R*"), above_1_Along = md("*C:R > 1*"), below_1_Along = md("*C:R < 1*")) %>%
    cols_hide("Function") %>%
    tab_footnote(footnote = "Includes parameter sets producing C:R > 10.",
        locations = cells_column_labels(columns = c("all_Against", "all_Along"))) %>%
    tab_footnote(footnote = "Does not include parameter sets producing C:R > 10.",
        locations = cells_column_labels(columns = c("above_1_Against",
            "above_1_Along"))) %>%
    tab_header(title = md("**Median LRR values of biomass and nutrient stocks**"),
        subtitle = md("LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing _C:R_ > 1, or only parameter sets producing _C:R_ < 1.")) %>%
    tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())

CRsummStock

Median LRR values of biomass and nutrient stocks
LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing C:R > 1, or only parameter sets producing C:R < 1.
Compartment	Against-gradient			Along-gradient
	All C:R1	C:R > 12	C:R < 1	All C:R1	C:R > 12	C:R < 1
Ecosystem 1
Consumers	−0.72	−0.71	−0.72	0.26	0.26	0.26
Nutrients	−0.16	−0.55	−0.15	0.12	0.20	0.08
Producers	0.00	0.00	0.00	0.00	0.00	0.00
Ecosystem 2
Consumers	0.17	0.15	0.10	−0.27	−0.13	−0.36
Nutrients	0.11	0.12	0.02	−0.14	−0.09	−0.10
Producers	0.04	0.03	0.10	−0.12	−0.12	−0.15
Meta-ecosystem
Consumers	0.08	−0.15	−0.12	−0.09	0.04	−0.09
Nutrients	0.03	−0.18	−0.04	−0.03	0.06	0.00
Producers	0.01	0.02	0.05	−0.02	−0.05	−0.04
1 Includes parameter sets producing C:R > 10.
2 Does not include parameter sets producing C:R > 10.

Nutrient Flux

Show code

CRsummFlux <- CRcomparison_tables %>%
    ungroup() %>%
    filter(., Function == "Nutrient Flux") %>%
    dplyr::mutate(., Compartment = fct_recode(Compartment, Consumers = "FLUX_C1rr",
        Consumers = "FLUX_C2rr", Consumers = "FLUX_Crr", Producers = "FLUX_P1rr",
        Producers = "FLUX_P2rr", Producers = "FLUX_Prr")) %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers",
        "Consumers")) %>%
    pivot_wider(., names_from = c(CR, Fertility), values_from = median) %>%
    gt(., groupname_col = "Scale") %>%
    fmt_number(., columns = 4:9, decimals = 2) %>%
    tab_spanner(label = md("**Along-gradient**"), columns = c("all_Along",
        "above_1_Along", "below_1_Along")) %>%
    tab_spanner(label = md("**Against-gradient**"), columns = c("all_Against",
        "above_1_Against", "below_1_Against")) %>%
    cols_label(Compartment = md("**Compartment**"), all_Against = md("*All C:R*"),
        above_1_Against = md("*C:R > 1*"), below_1_Against = md("*C:R < 1*"),
        all_Along = md("*All C:R*"), above_1_Along = md("*C:R > 1*"), below_1_Along = md("*C:R < 1*")) %>%
    cols_hide("Function") %>%
    tab_footnote(footnote = "Includes parameter sets producing C:R > 10.",
        locations = cells_column_labels(columns = c("all_Against", "all_Along"))) %>%
    tab_footnote(footnote = "Does not include parameter sets producing C:R > 10.",
        locations = cells_column_labels(columns = c("above_1_Against",
            "above_1_Along"))) %>%
    tab_header(title = md("**Median LRR values of nutrient flux**"), subtitle = md("LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing _C:R_ > 1, or only parameter sets producing _C:R_ < 1.")) %>%
    tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())

CRsummFlux

Median LRR values of nutrient flux
LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing C:R > 1, or only parameter sets producing C:R < 1.
Compartment	Against-gradient			Along-gradient
	All C:R1	C:R > 12	C:R < 1	All C:R1	C:R > 12	C:R < 1
Ecosystem 1
Consumers	−0.72	−0.71	−0.72	0.26	0.26	0.26
Producers	0.00	0.00	0.00	0.00	0.00	0.00
Ecosystem 2
Consumers	0.17	0.15	0.10	−0.27	−0.13	−0.36
Producers	0.04	0.03	0.10	−0.12	−0.12	−0.15
Meta-ecosystem
Consumers	0.07	−0.01	−0.07	−0.09	0.02	−0.11
Producers	0.01	0.02	0.04	−0.02	−0.05	−0.04
1 Includes parameter sets producing C:R > 10.
2 Does not include parameter sets producing C:R > 10.

Trophic Productivity

Show code

CRsummProd <- CRcomparison_tables %>%
    ungroup() %>%
    filter(., Function == "Trophic Productivity") %>%
    dplyr::mutate(., Compartment = fct_recode(Compartment, Consumers = "PROD_C1rr",
        Consumers = "PROD_C2rr", Consumers = "PROD_Crr", Producers = "PROD_P1rr",
        Producers = "PROD_P2rr", Producers = "PROD_Prr")) %>%
    dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers",
        "Consumers")) %>%
    pivot_wider(., names_from = c(CR, Fertility), values_from = median) %>%
    gt(., groupname_col = "Scale") %>%
    fmt_number(., columns = 4:9, decimals = 2) %>%
    tab_spanner(label = md("**Along-gradient**"), columns = c("all_Along",
        "above_1_Along", "below_1_Along")) %>%
    tab_spanner(label = md("**Against-gradient**"), columns = c("all_Against",
        "above_1_Against", "below_1_Against")) %>%
    cols_label(Compartment = md("**Compartment**"), all_Against = md("*All C:R*"),
        above_1_Against = md("*C:R > 1*"), below_1_Against = md("*C:R < 1*"),
        all_Along = md("*All C:R*"), above_1_Along = md("*C:R > 1*"), below_1_Along = md("*C:R < 1*")) %>%
    cols_hide("Function") %>%
    tab_footnote(footnote = "Includes parameter sets producing C:R > 10.",
        locations = cells_column_labels(columns = c("all_Against", "all_Along"))) %>%
    tab_footnote(footnote = "Does not include parameter sets producing C:R > 10.",
        locations = cells_column_labels(columns = c("above_1_Against",
            "above_1_Along"))) %>%
    tab_header(title = md("**Median LRR values of trophic compartment productivity**"),
        subtitle = md("LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing _C:R_ > 1, or only parameter sets producing _C:R_ < 1.")) %>%
    tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())

CRsummProd

Median LRR values of trophic compartment productivity
LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing C:R > 1, or only parameter sets producing C:R < 1.
Compartment	Against-gradient			Along-gradient
	All C:R1	C:R > 12	C:R < 1	All C:R1	C:R > 12	C:R < 1
Ecosystem 1
Consumers	−0.72	−0.71	−0.72	0.26	0.26	0.26
Producers	−0.16	−0.55	−0.15	0.12	0.20	0.08
Ecosystem 2
Consumers	0.22	0.20	0.22	−0.47	−0.28	−0.57
Producers	0.18	0.17	0.16	−0.34	−0.27	−0.34
Meta-ecosystem
Consumers	0.03	−0.07	−0.14	−0.04	0.06	−0.05
Producers	0.03	0.01	0.03	−0.04	−0.02	−0.05
1 Includes parameter sets producing C:R > 10.
2 Does not include parameter sets producing C:R > 10.

Sensitivity Analyses

Here, we perform a Global Sensitivity Analysis (GSA; sensu Harper, Stella, and Fremier (2011)) on the parameter sets that used to produce our results—i.e., those resulting in stable equilibria with biological meaning. Compared to local sensitivity analyses, whereby each parameter is varied over a range of values one at the time while all others are kept constant, GSA allows for varying all parameters at once. In turn, this has several advantages, namely (i) it removes the need to elect a constant value at which all parameters are kept while varying the one under analysis and (ii) it accounts for parameter interactions and is thus particularly suited for complex mathematical models with multiple parameters and non-linear relationships.

We use function randomForest() from package randomForest to run our sensitivity analysis. Rather than running a single analysis on the whole dataset, we run two separate ones—one each for the donor ecosystem 1 and the recipient ecosystem 2. This allows us to analyse parameter importance in each consumer movement scenario—gradient-neutral, along-, and against-gradient—and across the four different biomass distribution conditions identified above: (i) bottom-heavy pyramid in the donor ecosystem 1, (ii) bottom-heavy pyramid in the recipient ecosystem 2, (iii) top-heavy pyramid in the donor ecosystem 1, (iv) top-heavy pyramid in the recipient ecosystem 2

Setup

First, we load the ecoCR object, and create two additional variables: E1CRabove1 and E2CRabove1. Both these are categorical variables, with Yes or No values, and capture whether the local biomass distribution in the donor ecosystem 1 or in the recipient ecosystem 2 conforms to a classic bottom-heavy pyramid or to an inverted, top-heavy pyramid.

Afterwards, we store the names of the six state variables in our model into a new object stateV that we will use later in the code to isolate each state variable in turn.

Show code

# load the data, then create additional categorical variables for
# plotting
ecoCR <- readRDS(file = "../Results/ecoCR.rds") %>%
    mutate(., E1CRabove1 = as_factor(ifelse(Eco1 > 1, "Y", "N")), E2CRabove1 = as_factor(ifelse(Eco2 >
        1, "Y", "N")), Fertility = as_factor(Fertility))

# define a vector to use in the for loop to go through each state
# variable in the model
stateV <- c("N1", "P1", "C1", "N2", "P2", "C2")

Sensitivity Analyses by Spatial Scale

We will run two Global Sensitivity Analyses on the whole ecoCR object, based on whether we take into account the biomass distribution in the donor ecosystem 1 or in the recipient ecosystem 2. This is necessary as the same parameter set can give rise to an inverted, top-heavy biomass pyramid in either local ecosystem but to a classic, bottom-heavy one in the other. By focusing on one ecosystem at the time in terms of biomass distribution, we can account for all four scenarios highlighted above.

Donor ecosystem 1

Starting with the donor ecosystem 1, we use a for loop to run out GSA on our ecoCR object, one state variable at the time. For each run of the loop, we isolate the relevant state variable, generate the random forest trees, calculate importance and relative importance, and store all this in a new object, sensiResE1.

Results for this analysis are shown in Figures 27 and 29.

Show code

# allocate an empty dataframe to contain the results of each
# sensitivity analysis
sensiResE1 <- NULL

# This loop runs the sensitivity analysis
for (i in 1:length(stateV)) {

    # wrangle the dataframe to isolate the columns we want to work
    # with Note that we are only include E1CRabove1 because here we
    # focus on biomass distribution (i.e., C:R value) in the donor
    # ecosystem 1
    sensiA <- select(ecoCR, stateV[i], I1:e2, Fertility, E1CRabove1)

    # rename the state variable column to a neutral name
    colnames(sensiA)[1] <- "Compartment"

    # here I group by Fertility (= movement scenario) and whether the
    # CR in Ecosystem 1 (donor) is > 1, so I can later use these two
    # pieces of information to produce multiple graphs from the same
    # df (which reduces the number of runs for the for loop) The last
    # function call in the pine runs the randomForest model to
    # estimate
    sA_temp <- sensiA %>%
        group_by(Fertility, E1CRabove1) %>%
        select(., Compartment:e2) %>%
        do(Cmove = randomForest(Compartment ~ . - Fertility - E1CRabove1,
            data = ., proximity = TRUE, importance = TRUE))  # run this and see

    # allocate empty dataframe to store results from next for loop
    imp_temp <- NULL

    # this loop gets the results for each movement scenario
    for (j in seq(1, 6, 1)) {
        # this pipe extracts the importance values of each parameter
        # and stores them in the df created just above
        temp <- as_tibble(importance(sA_temp[[3]][[j]]), rownames = NA) %>%
            rownames_to_column(., var = "Parameter") %>%
            mutate(., Movement = as_factor(sA_temp[[1]][[j]]), E1CRabove1 = as_factor(sA_temp[[2]][[j]]))
        imp_temp <- rbind(imp_temp, temp)
        # browser() print(j)
    }

    # now calculate the ranked importance value as the ratio of each
    # parameter's importance values to the sum of all parameters'
    # importance values, store it in a temporary object to which we
    # add grouping and categorical variables that are useful in
    # plotting the results
    ranked_imp_temp <- tibble(imp_temp[, 3]/sum(imp_temp[, 3])) %>%
        mutate(., Parameter = as_factor(imp_temp$Parameter), Movement = as_factor(imp_temp$Movement),
            E1CRabove1 = as_factor(imp_temp$E1CRabove1), SV = as_factor(stateV[i]))

    # store the results from each loop run in a separate df so the
    # loop can run again without loosing results
    sensiResE1 <- bind_rows(sensiResE1, ranked_imp_temp)

    # browser()

    print(i)

}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6

Show code

# remove objects used in the loop to free up memory and avoid errors
# due to values being transferred over to the next sensitivity
# analysis
rm(list = c("i", "j", "imp_temp", "sA_temp", "ranked_imp_temp", "temp",
    "sensiA"))

Recipient ecosystem 2

The GSA taking into account the biomass distribution of ecosystem 2 runs along the same lines as the one for the donor ecosystem 1, with the results being stored in object sensiResE2.

Results for this analysis are shown in Figures 28 and 30.

Show code

# allocate an empty dataframe to contain the results of each
# sensitivity analysis
sensiResE2 <- NULL

# This loop runs the sensitivity analysis
for (i in 1:length(stateV)) {

    # wrangle the dataframe to isolate the columns we want to work
    # with Note that we are only include E2CRabove1 because here we
    # focus on biomass distribution (i.e., C:R value) in the
    # recipient ecosystem 2
    sensiA <- select(ecoCR, stateV[i], I1:e2, Fertility, E2CRabove1)

    # rename the state variable column to a neutral name
    colnames(sensiA)[1] <- "Compartment"

    # here I group by Fertility (= movement scenario) and whether the
    # CR in Ecosystem 2 (recipient) is > 1, so I can later use these
    # two pieces of information to produce multiple graphs from the
    # same df (which reduces the number of runs for the for loop) The
    # last function call in the pine runs the randomForest model to
    # estimate
    sA_temp <- sensiA %>%
        group_by(Fertility, E2CRabove1) %>%
        select(., Compartment:e2) %>%
        do(Cmove = randomForest(Compartment ~ . - Fertility - E2CRabove1,
            data = ., proximity = TRUE, importance = TRUE))  # run this and see

    # allocate empty dataframe to store results from next for loop
    imp_temp <- NULL

    # this loop gets the results for each movement scenario
    for (j in seq(1, 6, 1)) {
        # this pipe extracts the importance values of each parameter
        # and stores them in the df created just above
        temp <- as_tibble(importance(sA_temp[[3]][[j]]), rownames = NA) %>%
            rownames_to_column(., var = "Parameter") %>%
            mutate(., Movement = as_factor(sA_temp[[1]][[j]]), E2CRabove1 = as_factor(sA_temp[[2]][[j]]))
        imp_temp <- rbind(imp_temp, temp)
        # browser() print(j)
    }

    # now calculate the ranked importance value as the ratio of each
    # parameter's importance values to the sum of all parameters'
    # importance values, store it in a temporary object to which we
    # add grouping and categorical variables that are useful in
    # plotting the results
    ranked_imp_temp <- tibble(imp_temp[, 3]/sum(imp_temp[, 3])) %>%
        mutate(., Parameter = as_factor(imp_temp$Parameter), Movement = as_factor(imp_temp$Movement),
            E2CRabove1 = as_factor(imp_temp$E2CRabove1), SV = as_factor(stateV[i]))

    # store the results from each loop run in a separate df so the
    # loop can run again without loosing results
    sensiResE2 <- bind_rows(sensiResE2, ranked_imp_temp)

    # browser()

    print(i)

}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6

Show code

# remove objects used in the loop to free up memory and avoid errors
# due to values being transferred over to the next sensitivity
# analysis
rm(list = c("i", "j", "imp_temp", "sA_temp", "ranked_imp_temp", "temp",
    "sensiA"))

Whole system

As shown in the figures below (Figures @ref(fig: gsa-results-CR-below-1-eco-1) to 30), the shape of the biomass pyramid in either the donor or recipient ecosystem does not appear to influence the relative importance of the parameters in the model. For this reason, we run a final GSA on the whole model. Results from this model run are shown in Figure 26.

Show code

# allocate an empty dataframe to contain the results of each
# sensitivity analysis
sensiRes <- NULL

# This loop runs the sensitivity analysis
for (i in 1:length(stateV)) {

    # wrangle the dataframe to isolate the columns we want to work
    # with Note that we are only include E2CRabove1 because here we
    # focus on biomass distribution (i.e., C:R value) in the
    # recipient ecosystem 2
    sensiA <- select(ecoCR, stateV[i], I1:e2, Fertility)

    # rename the state variable column to a neutral name
    colnames(sensiA)[1] <- "Compartment"

    # here I group by Fertility (= movement scenario) and whether the
    # CR in Ecosystem 2 (recipient) is > 1, so I can later use these
    # two pieces of information to produce multiple graphs from the
    # same df (which reduces the number of runs for the for loop) The
    # last function call in the pine runs the randomForest model to
    # estimate
    sA_temp <- sensiA %>%
        group_by(Fertility) %>%
        select(., Compartment:e2) %>%
        do(Cmove = randomForest(Compartment ~ . - Fertility, data = .,
            proximity = TRUE, importance = TRUE))  # run this and see

    # allocate empty dataframe to store results from next for loop
    imp_temp <- NULL

    # this loop gets the results for each movement scenario
    for (j in seq(1, 3, 1)) {
        # this pipe extracts the importance values of each parameter
        # and stores them in the df created just above
        temp <- as_tibble(importance(sA_temp[[2]][[j]]), rownames = NA) %>%
            rownames_to_column(., var = "Parameter") %>%
            mutate(., Movement = as_factor(sA_temp[[1]][[j]]))
        imp_temp <- rbind(imp_temp, temp)
        # browser() print(j)
    }

    # now calculate the ranked importance value as the ratio of each
    # parameter's importance values to the sum of all parameters'
    # importance values, store it in a temporary object to which we
    # add grouping and categorical variables that are useful in
    # plotting the results
    ranked_imp_temp <- tibble(imp_temp[, 3]/sum(imp_temp[, 3])) %>%
        mutate(., Parameter = as_factor(imp_temp$Parameter), Movement = as_factor(imp_temp$Movement),
            SV = as_factor(stateV[i]))

    # store the results from each loop run in a separate df so the
    # loop can run again without loosing results
    sensiRes <- bind_rows(sensiRes, ranked_imp_temp)

    # browser()

    print(i)

}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6

Show code

# remove objects used in the loop to free up memory and avoid errors
# due to values being transferred over to the next sensitivity
# analysis
rm(list = c("i", "j", "imp_temp", "sA_temp", "ranked_imp_temp", "temp",
    "sensiA"))

Visual Deliverables—Relative Parameter Importance

Now that we have a relative importance value for each parameter in each movement scenario, in all four conditions of biomass distribution in the two local ecosystems (see above, Changes to biomass distribution in local and meta-ecosystem), we produce graphs to visualize these results. These graphs will be useful in identifying the more important parameters in determining the stock value of each trophic compartment in the model for each combination of movement scenario and biomass distribution conditions.

Figure 26 presents the results from the GSA run on the whole system, which do not differ from those of two GSA runs on the two local ecosystems (see above; Whole System). Figures 27 to 30 below present the results from the local-ecosystem GSAs, for transparancy and completeness.

Whole system

Show code

# plotting pipe first, rename the variable to use as x-axis, then
# filter by C:R > 1 in the recipient ecosystem 2, then create a new
# variable to store information about which ecosystem do parameters
# belong to. Then, turn it into a factor, reorder the levels using
# fct_relevel(), and reorder the levels of factor Parameter and their
# associated relative importance values according to the Scale
# factor.  Finally, group by Scale and plot
GSAwhole <- sensiRes %>%
    rename(., RelImpI = IncNodePurity) %>%
    mutate(., Scale = ifelse(Parameter == "u1" | Parameter == "h1" | Parameter ==
        "a1" | Parameter == "e1" | Parameter == "d1" | Parameter == "I1",
        "Donor", ifelse(Parameter == "g" | Parameter == "m" | Parameter ==
            "c", "Dispersers' Pool", ifelse(Parameter == "l", "Other",
            "Recipient"))), Scale = as_factor(Scale), Scale = fct_relevel(Scale,
        "Donor", "Recipient", "Dispersers' Pool", "Other"), Parameter = fct_reorder2(.f = Parameter,
        .x = RelImpI, .y = Scale)) %>%
    group_by(., Scale) %>%
    ggplot(., aes(x = RelImpI, y = Parameter)) + geom_col(aes(fill = Scale),
    col = "black", linewidth = 0.25) + scale_fill_met_d(name = "Egypt",
    direction = 1) + xlab("Relative Importance Index") + ylab(" ") + facet_grid2(Movement ~
    SV) + theme(legend.position = "top", legend.title = element_blank(),
    text = element_text(size = 20), panel.spacing = unit(0.5, "lines"))

GSAwhole + scale_y_discrete("Parameter", labels = c("l", "g", "m", "c",
    expression(I[2]), expression(u[2]), expression(a[2]), expression(h[2]),
    expression(d[2]), expression(e[2]), expression(I[2]), expression(u[1]),
    expression(a[1]), expression(h[1]), expression(d[1]), expression(e[1])))

Figure 26: At the whole system level, considering together all cases of biomass distribution in the two ecosystems, the global sensitivity analysis again shows evidence for the direct and indirect effects of functional response-related parameters across all movement scenarios and in both ecosystems. In particular, when consumers move against-gradient (bottom row), notice the importance of parameters e2 and d2 to shape not only the dynamics of trophic compartment in the recipient ecosystem 2 but also those of the primary producers in ecosystem 1. Note that this figure is equivalent to Figure A.7 in the Appendices to the main text, but with x-axes tailored to each panel’s range of relative importance values.

Show code

# save statement
ggsave("../Results/SA_RelImpI_WS.pdf", device = "pdf", dpi = 300, width = 12,
    height = 10)

Classic donor

Show code

# plotting pipe first, rename the variable to use as x-axis, then
# filter by C:R < 1 in the donor ecosystem 1, then create a new
# variable to store information about which ecosystem do parameters
# belong to. Then, turn it into a factor, reorder the levels using
# fct_relevel(), and reorder the levels of factor Parameter and their
# associated relative importance values according to the Scale
# factor.  Finally, group by Scale and plot
sensiResE1 %>%
    rename(., RelImpI = IncNodePurity) %>%
    filter(., E1CRabove1 == "N") %>%
    mutate(., Scale = ifelse(Parameter == "u1" | Parameter == "h1" | Parameter ==
        "a1" | Parameter == "e1" | Parameter == "d1" | Parameter == "I1",
        "Donor", ifelse(Parameter == "g" | Parameter == "m" | Parameter ==
            "c", "Dispersers' Pool", ifelse(Parameter == "l", "Other",
            "Recipient"))), Scale = as_factor(Scale), Scale = fct_relevel(Scale,
        "Donor", "Recipient", "Dispersers' Pool", "Other"), Parameter = fct_reorder2(.f = Parameter,
        .x = RelImpI, .y = Scale)) %>%
    group_by(., Scale) %>%
    ggplot(., aes(x = RelImpI, y = Parameter)) + geom_col(aes(fill = Scale),
    col = "black", linewidth = 0.25) + scale_fill_met_d(name = "Egypt",
    direction = 1) + xlab("Relative Importance Index") + ylab(" ") + facet_grid2(Movement ~
    SV) + theme(legend.position = "top", legend.title = element_blank(),
    text = element_text(size = 20), panel.spacing = unit(1.5, "lines")) +
    scale_y_discrete("Parameter", labels = c("l", "g", "m", "c", expression(I[2]),
        expression(u[2]), expression(a[2]), expression(h[2]), expression(d[2]),
        expression(e[2]), expression(I[2]), expression(u[1]), expression(a[1]),
        expression(h[1]), expression(d[1]), expression(e[1])))

Figure 27: The global sensitivity analysis (GSA) performed on the subset of parameter sets that produces bottom-heavy biomass pyramids in the donor ecosystem provides key insights in the direct and indirect pathways through which higher trophic compartment influence lower ones. The high relative importance of parameter u1, involved in the functional response of Primary Producers, is a clear example of a direct effect shaping the accumulation of Nutrients in the donor ecosystem 1 in all three movement scenario (first column). Likewise, the high relative importance of parameter e2, the Consumer efficiency in the recipient ecosystem 2, in shaping Nutrient stocks in the recipient ecosystem 2 (fourth column, bottom row) highlights the strength of the indirect pathway that allows the highest trophic level in our system to influence the lower one. Note that e2 is also the most importance parameter in determining biomass accumulation of Consumers in the recipient ecosystem 2 in this scenario, meaning that it has both direct and indirect effects on the system’s dynamics.

Show code

# save statement
ggsave("../Results/SA_RelImpI_Eco1CRbelow1.pdf", device = "pdf", dpi = 300,
    width = 12, height = 10)

Classic recipient

Show code

# plotting pipe first, rename the variable to use as x-axis, then
# filter by C:R < 1 in the recipient ecosystem 2, then create a new
# variable to store information about which ecosystem do parameters
# belong to. Then, turn it into a factor, reorder the levels using
# fct_relevel(), and reorder the levels of factor Parameter and their
# associated relative importance values according to the Scale
# factor.  Finally, group by Scale and plot
sensiResE2 %>%
    rename(., RelImpI = IncNodePurity) %>%
    filter(., E2CRabove1 == "N") %>%
    mutate(., Scale = ifelse(Parameter == "u1" | Parameter == "h1" | Parameter ==
        "a1" | Parameter == "e1" | Parameter == "d1" | Parameter == "I1",
        "Donor", ifelse(Parameter == "g" | Parameter == "m" | Parameter ==
            "c", "Dispersers' Pool", ifelse(Parameter == "l", "Other",
            "Recipient"))), Scale = as_factor(Scale), Scale = fct_relevel(Scale,
        "Donor", "Recipient", "Dispersers' Pool", "Other"), Parameter = fct_reorder2(.f = Parameter,
        .x = RelImpI, .y = Scale)) %>%
    group_by(., Scale) %>%
    ggplot(., aes(x = RelImpI, y = Parameter)) + geom_col(aes(fill = Scale),
    col = "black", linewidth = 0.25) + scale_fill_met_d(name = "Egypt",
    direction = 1) + xlab("Relative Importance Index") + ylab(" ") + facet_grid2(Movement ~
    SV) + theme(legend.position = "top", legend.title = element_blank(),
    text = element_text(size = 20), panel.spacing = unit(1.5, "lines")) +
    scale_y_discrete("Parameter", labels = c("l", "g", "m", "c", expression(I[2]),
        expression(u[2]), expression(a[2]), expression(h[2]), expression(d[2]),
        expression(e[2]), expression(I[2]), expression(u[1]), expression(a[1]),
        expression(h[1]), expression(d[1]), expression(e[1])))

Figure 28: When the recipient ecosystem 2 has a bottom-heavy biomass pyramid, we see evidence of direct and indirect effects of high trophic compartments on lower ones arising from our global sensitivity analysis. Of particular note are the direct and indirect effects of parameter e1 on the biomass accumulation of Consumers and Nutrient in the donor ecosystem 1, respectively, when consumer move along-gradient.

Show code

# save statement
ggsave("../Results/SA_RelImpI_Eco2CRbelow1.pdf", device = "pdf", dpi = 300,
    width = 12, height = 10)

Inverted donor

Show code

# plotting pipe first, rename the variable to use as x-axis, then
# filter by C:R > 1 in the donor ecosystem 1, then create a new
# variable to store information about which ecosystem do parameters
# belong to. Then, turn it into a factor, reorder the levels using
# fct_relevel(), and reorder the levels of factor Parameter and their
# associated relative importance values according to the Scale
# factor.  Finally, group by Scale and plot
sensiResE1 %>%
    rename(., RelImpI = IncNodePurity) %>%
    filter(., E1CRabove1 == "Y") %>%
    mutate(., Scale = ifelse(Parameter == "u1" | Parameter == "h1" | Parameter ==
        "a1" | Parameter == "e1" | Parameter == "d1" | Parameter == "I1",
        "Donor", ifelse(Parameter == "g" | Parameter == "m" | Parameter ==
            "c", "Dispersers' Pool", ifelse(Parameter == "l", "Other",
            "Recipient"))), Scale = as_factor(Scale), Scale = fct_relevel(Scale,
        "Donor", "Recipient", "Dispersers' Pool", "Other"), Parameter = fct_reorder2(.f = Parameter,
        .x = RelImpI, .y = Scale)) %>%
    group_by(., Scale) %>%
    ggplot(., aes(x = RelImpI, y = Parameter)) + geom_col(aes(fill = Scale),
    col = "black", linewidth = 0.25) + scale_fill_met_d(name = "Egypt",
    direction = 1) + xlab("Relative Importance Index") + ylab(" ") + facet_grid2(Movement ~
    SV) + theme(legend.position = "top", legend.title = element_blank(),
    text = element_text(size = 20), panel.spacing = unit(1.5, "lines")) +
    scale_y_discrete("Parameter", labels = c("l", "g", "m", "c", expression(I[2]),
        expression(u[2]), expression(a[2]), expression(h[2]), expression(d[2]),
        expression(e[2]), expression(I[2]), expression(u[1]), expression(a[1]),
        expression(h[1]), expression(d[1]), expression(e[1])))

Figure 29: When the donor ecosystem 1 has a top-heavy biomass pyramid, we note again strong direct effects of functional response-related parameters on the dynamics of Consumers and Nutrient in this ecosystem, when consumers move either in a gradient-neutral or along-gradient fashion. Note that weak indirect effects of parameter e1 on Nutrient stock accumulation in the donor ecosystem 1 are still visible when consumers move along-gradient (first column, middle row). In the against-gradient movement scenario, we note again some evidence of indirect effects of parameter e2 on the dynamics of the Nutrient compartment in the recipient ecosystem 2 (fourth column, bottom row).

Show code

# save statement
ggsave("../Results/SA_RelImpI_Eco1CRabove1.pdf", device = "pdf", dpi = 300,
    width = 12, height = 10)

Inverted recipient

Show code

# plotting pipe first, rename the variable to use as x-axis, then
# filter by C:R > 1 in the recipient ecosystem 2, then create a new
# variable to store information about which ecosystem do parameters
# belong to. Then, turn it into a factor, reorder the levels using
# fct_relevel(), and reorder the levels of factor Parameter and their
# associated relative importance values according to the Scale
# factor.  Finally, group by Scale and plot
sensiResE2 %>%
    rename(., RelImpI = IncNodePurity) %>%
    filter(., E2CRabove1 == "Y") %>%
    mutate(., Scale = ifelse(Parameter == "u1" | Parameter == "h1" | Parameter ==
        "a1" | Parameter == "e1" | Parameter == "d1" | Parameter == "I1",
        "Donor", ifelse(Parameter == "g" | Parameter == "m" | Parameter ==
            "c", "Dispersers' Pool", ifelse(Parameter == "l", "Other",
            "Recipient"))), Scale = as_factor(Scale), Scale = fct_relevel(Scale,
        "Donor", "Recipient", "Dispersers' Pool", "Other"), Parameter = fct_reorder2(.f = Parameter,
        .x = RelImpI, .y = Scale)) %>%
    group_by(., Scale) %>%
    ggplot(., aes(x = RelImpI, y = Parameter)) + geom_col(aes(fill = Scale),
    col = "black", linewidth = 0.25) + scale_fill_met_d(name = "Egypt",
    direction = 1) + xlab("Relative Importance Index") + ylab(" ") + facet_grid2(Movement ~
    SV) + theme(legend.position = "top", legend.title = element_blank(),
    text = element_text(size = 20), panel.spacing = unit(1.5, "lines")) +
    scale_y_discrete("Parameter", labels = c("l", "g", "m", "c", expression(I[2]),
        expression(u[2]), expression(a[2]), expression(h[2]), expression(d[2]),
        expression(e[2]), expression(I[2]), expression(u[1]), expression(a[1]),
        expression(h[1]), expression(d[1]), expression(e[1])))

Figure 30: When the recipient ecosystem 2 has a top-heavy biomass pyramid, the high relative importance of functional response-related parameters in both the dynamics of the Consumer and Nutrient compartments of the recipient ecosystem 2 points to both direct and indirect effects being exerted by Consumers in the against-gradient and gradient-neutral movement scenarios (forth and sixth columns, top and bottom rows).

Show code

# save statement
ggsave("../Results/SA_RelImpI_Eco2CRabove1.pdf", device = "pdf", dpi = 300,
    width = 12, height = 10)

Compilation Time

This notebook was compiled in 37 minutes.

Acknowledgments

We are grateful to Dr. Anne McLeod for her help in developing the stability and sensitivity analyses code, and for suggesting useful reference on both topics. We thank Dr. Robert Buchkowski for his insights in designing and evaluating the stability analyses. Finally, we thank Dr. A. Z. Andis Arietta for his help in developing the code for the biomass pyramid plots.

Author Contributions

MR and SJL conceptualized the model; MR, SJL, OJS, EVW, YFW, and TRH designed the study; MR wrote the code and run the analyses; MR and SJL run and proofread the code; MR, SJL, OJS interpreted the results; all authors contributed critically to the manuscript.